ARLTR-85-32 


Copy  No. 


NONLINEAR  EFFECTS  IN  LONG  RANGE 
UNDERWATER  ACOUSTIC  PROPAGATION 


Frederick  0.  Coteras 


APPLIED  RESEARCH  LABORATORIES 

THE  UNIVERSITY  OF  TEXAS  AT  AUSTIN 
FOST  OFFICE  BOX  8029,  AUSTIN.  TEXAS  78713-8029 


1  November  198S 
Technical  Report 


Approved  for  public  release; 
distribution  unlimited. 


OTIC 

tVLECTE 

APR08«6 


Prepared  for: 


OFFICE  OF  NAVAL  RESEARCH 
DEPARTMENT  OF  THE  NAVY 
ARLINGTON,  VA  22217 


/ 


UNCLASSIFIED 

security  classification  of  this  pace  rllflwi  Data  eniatad) 


REPORT  DOCUMENTATION  PAGE 


REPORT  number 


4.  title  (mN  SubHtIa) 


TT’nrmT'Tra 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


RECIPIENT'S  CATALOG  NUMBER 


NONLINEAR  EFFECTS  IN  LONG  RANGE  UNDERWATER 
ACOUSTIC  PROPAGATION 


7.  author^*; 


Frederick  D.  Cotaras 


9.  PERrORMING  ORGANIZATION  NAME  AND  AOORESS 

Applied  Research  Laboratories 
The  University  of  Texas  at  Austin 
Austin,  Texas  78713-8029 


n.  CONTROLLING  OFFICE  NAME  AND  AOORESS 


s.  Type  of  report  a  perioo  covered 

technical  report 


s.  performing  org.  report  number 

ARL-TR-85-32 


S.  CONTRACT  OR  GRANT  NUMBERf*; 

N00014-75-C-0867 

N00014-84-K-0574 


Office  of  Naval  Research 
Department  of  the  Navy 
Arlington,  Virginia  22217 


12.  report  DATE 

1  November  1985 


13.  NUM0ER  OF  PAGES 

217 


MONITORING  AGENCY  NAME  B  AOORESSfli  dtUatant  /ran  Cantralllnt  OUlea)  IS.  SECURITY  CLASS,  (al  thim  nport) 

UNCLASSIFIED 

ISa.^ECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (al  iMa  Rapotl) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (al  lha  abalrael  aniarad  In  BlaeM  20,  II  dlllarani  from  Report; 


19.  KEY  WORDS  (Conllnaa  an  rararaa  alda  II  nacaaaary  and  Idantlly  by  black  nuitibar) 


underwater  explosions  numerical  model 

ray  theory  finite  amplitude  effects 

long  range  underwater  propagation 
eikonal  equation 


20.  ABSTRACT  (Conllnua  on  raaaraa  alda  II  nacaaaary  and  Idantlly  by  black  numbar) 

In  this  report '^tbe  propagation  of  finite  amplitude  acoustic  signals  through  an 
inhomogeneous  ocean  ‘fs- invest iga^t€^ both  analytically  and  numerically.  The 
effects  of  reflections  and  focusing  are  not  considered.  From  simplified 
versions  of  the  lossless  hydrodynamics  equations  the  theories  of  linear  and 
nonlinear  geometrical  acoustics  are  developed.  Losses  are  accounted  for 
directly  in  the  numerical  routine.  The  eikonal  equation,  from  which  an  equation 
for  the  ray  paths  is  derived,  is  assumed  to  be  the  same  for  both  small -signal 


I  JAN  73  1473  eOlTlON  OF  1  NOV  66  IS  OBSOLETE  UNCLASSIFIED 

f  A  security  classification  of  This  PAGEfBNwi  Dm*  BntpraW; 


UNCLASSIFIED 


SKCUWTY  CUMUFICATIOM  or  TMII  Of  Kmifd) 


20.  (cont'd) 

and  finite  amplitude  waves.  The  transport  equation  is  found  to  be  different, 
however.  The  transport  equation  leads  to  a  standard  first-order  progressive 
wave  equation,  linear  for  small-signals  waves,  but  nonlinear  for  finite 
amplitude  waves.  All  the  analysis  is  carried  out  in  the  time  domain  and  is 
for  a  fully  Inhomogeneous  ocean. 

>In  the  numerical  study  the  ocean  is  assumed  to  be  stratified.  The  effects  of 
inhomogeneity,  ordinary  attenuation  and  dispersion,  and  nonlinear  propagation 
are  investigated  using  a  numerical  implementation  of  nonlinear  geometrical 
acoustics.  Two  explosion  waveforms  are  considered;  a  weak  shock  with  an 
exponentially  decaying  tail  and  a  more  realistic  waveform  that  includes  the 
first  bubble  pulse.  Numerical  propagation  of  the  simpler  wave  along  a  58.1  km 
path  starting  at  a  depth  of  300  m  leads  to  the  following  conclusions:  (1)  The 
effept  of  inhomogeneity  on  nonlinear  distortion  is  small.  (2)  Dispersion  plays 
anr  important  role  in  determining  the  arrival  time  of  the  pulse.  (3)  Neither 
nonlinearity  nor  ordinary  attenuation  (and  dispersion)  are  paramount;  both 
need  tp  be  included.  For  the  more  realistic  wave propagation  is  along  a 
23  km  ray  path  starting  from  a  depth  of  4300  m.  Two  charge  weights,  0.818  kg 
and  22.7  kg  TNT,  are  assumed.  In  each  case  the  energy  spectrum  of  the  signal 
obtained  by  considering  finite  amplitude  effects  for  the  entire  23  km  path 
is  compared  with  spectra  obtained  by  neglecting  finite  amplitude  effects 
(li  entirely,  42)  after  the  first  150  m,  and  4^ after  the  first  1100  m. 

Finite  amplitude  effects  are^fb«n€Tto  b^-of  small  consequence  in  the  case  of 
the  0.818  kg  TNT  explosion  for  frequencies  below  6  kHz  at  distances  beyond 
1100  m.  For  the  22.7  kg  explosion  the  corresponding  quantities  are  4  kHz 
and  1100  m.  ^  f/c-'"  ^.<y  i  "i  ■'■<  (  i  ■ 

I  V  ■' 

)  >  j'  .  .  ^  I  '• '  • 

t 


UNCLASSIFIED 

StCURITV  CkAtSirtCATlON  OP  THIS  RAOKVDmi  Data  iMtara^ 


TAEU.E  OF  CONTEHTS 


Eaflfi 


LIST  OF  FIGURES  5 

LIST  OF  IMPORTANT  SYMBOLS 

FOREWORD  13 

CHAPTER  1  INTRODUaiON  15 

A.  Review  of  Finite  Amplitude  Effects  1 7 

B.  Review  of  Linear  Geometrical  Acoustics  19 

C.  Geometrical  Acoustics  for  Finite  Amplitude  Signals  23 

0.  Scope  of  the  Investigation  26 

CHAPTER  2  RANKING  OF  TERMS  28 

A.  Introduction  28 

B.  Simplifying  the  Lossless  Hydrodynamics  Equations  33 

1.  Linear  Hydrodynamics  Equations  for  Lossless 

Homogeneous  Fluids  34 

2.  Linear  Hydrodynamics  Equations  for  Lossless 

Inhomogeneous  Fluids  35 

T  Accesion  For  \  I 


3.  Nonlinear  Hydrodynamics  Equations  for  Lossless 
Inhomogeneous  Fluids 

C,  Substitution  into  Second-Order  Terms  Using  First-Order 
Relations 

D.  Epilogue 

CHAPTERS  LINEAR  GEOMETRICAL  ACOUSTICS 
A  Introduction 

B.  Linear  Lossless  Wave  Equation  for  Inhomogeneous  Fluids 

C.  Geometrical  Acoustics  Assumption 

D.  Ray  Paths  from  the  Eikonal 

E.  Acoustic  Pressure  from  the  Transport  Equation 
CHAPTER  4  NONLINEAR  GEOrtTRICAL  ACOUSTICS 

A  Combining  the  Hydrodynamics  Equations 

B.  Geometrical  Acoustics  Assumption 

C  Reduction  of  the  Nonlinear  Geometrical  Acoustics 
Equation 

D.  The  Solution 

E.  Simplification  of  the  Distortion  Range  Variable 
Transformation 

CHAPTERS  NUMERICAL  EVALUATION  OF  WAVEFORMS 
A  Introduction 

B.  Application  of  Pestorius's  Algorithm  to  Inhomogeneous 
Media 


I.  The  Weak  Shock  Subroutine 


3 

2m 


2.  The  Appllcatton  of  Attenuation  84 

3.  Other  Considerations  90 

C.  Verification  of  the  Algorithm  92 

CHAPTER  6  RESULTS  100 

A  Introduction  100 

B.  Design  of  the  Numerical  Experiment  101 

1.  Ray  Paths  101 

2.  Signals  106 

C.  Effect  of  Medium  Inhomogeneity  on  Nonlinearity  1 10 

D.  The  Combined  Effects  of  Nonlinearity  and  Absorption  1 1 4 


E.  To  What  Distance  are  Finite  Amplitude  Effects  Important?  1 20 

1 .  The  0.81 8  kg  TNT  Explosion  Results 

2.  The  22.7  kg  TNT  Explosion  Results 
CHAPTER  7  SUMMARY 

APPENDIX  A  RAY  COORDINATES  AND  THE  EXPRESSION  FOR 

APPENDIX  B  THE  PARAMETER  OF  NONLINEARITY  FOR  SEAWATER 

APPENDIX  C  ANALYTIC  SaUTIONS  FOR  FINITE  AMPLITUDE  WAVES 
VIA  WEAK  SHOCK  THEORY 

APPENDIX  D  COMPUTER  PROGRAM 

REFERENCES 


121 

126 

132 

136 

144 

153 

163 

201 


LIST  OF  FI6URES 


Figure 

IlUfi 

eag£ 

3.1 

Ray  Paths  and  Equiphase  Wavefronts 

43 

3.2 

The  Ray  Coordinate  System 

52 

5.1 

Flowchart  of  the  Modified  Pestorlus  Algorithm 

81 

5.2 

Salinity  Profile  Typical  of  North  Atlantic  Ocean 

85 

5.3 

Temperature  Profile  Typical  of  North  Atlantic  Ocean 

86 

5.4 

Time  Waveforms  for  Exponentially  Decaying  Pulse 
Propagating  Spherically  through  a  Lossless  Ocean 

97 

5.5 

Energy  Spectra  for  Exponentially  Decaying  Pulse  Propagating 
Spherically  through  a  Lossless  Homogeneous  Ocean 

98 

6.1 

Ocean  Environment  and  Acoustic  Ray  Paths 

102 

6.2 

Family  of  Rays  from  0^  to  15*  at  0.5®  Increments  Starting 
from  300  m 

104 

6.3 

Family  of  Rays  from  5®  to  10®  at  0.5®  Increments  Starting 
from  4300  m 

105 

6.4 

Time  Waveforms  and  Energy  Spectra  for  300  m  Detonation 

109 

6.5 

Time  Waveforms  and  Energy  Spectra  for  4300  m  Detonation 

ill 

Weak  Shock  with  an  Exponentially  Decaying  Tail  after 
Propagating  58. 1  km  through  ( 1 )  Lossless  Stratified 
Ocean  and  (2)  Lossless  Homogeneous  Ocean 

Weak  Shock  with  an  Exponentially  Decaying  Tail  after 
Propagating  58. 1  km  through  a  Stratified  Ocean 
Considering  ( 1 )  Ordinary  Attenuation  and  Dispersion  only 
(2)  Finite  Amplitude  Effects  only  and  (3)  Finite  Amplitude 
Effects  and  Ordinary  Attenuation  and  Dispersion 

Propagation  of  a  Weak  Shock  with  an  Exponentially 
Decaying  Tail  through  a  Lossy  Stratified  Ocean 
Neglecting  Finite  Amplitude  Effects 

Propagation  of  a  Weak  Shock  with  an  Exponentially 
Decaying  Tail  through  a  Lossy  Stratified  Ocean 
Considering  Finite  Amplitude  Effects 

Energy  Spectra  of  0.8 1 8  kg  TNT  Explosion  ( I )  at  the 
Reference  Range  (0.4  m),  (2)  after  Propagating  23  km 
Neglecting  Finite  Amplitude  Effects  (Case  A),  and  (3) 
after  Propagating  23  km  Considering  Finite  Amplitude 
Effects  (Case  D) 

Energy  Spectra  of  0.818  kg  TNT  Explosion  after  Propagating 
23  km  Considering  Finite  Amplitude  Effects  ( I )  over  the 
First  150  m  (Case  8)  and  (2)  over  the  Entire  23  km  (Case  D) 

Energy  Spectra  of  0.818  kg  TNT  Explosion  after  Propagating 
23  km  Considering  Finite  Amplitude  Effects  over  ( I ) 
the  First  1  iOO  m  (Case  C)  and  (2)  the  Entire  Propagation 
Path  (Case  D) 

Energy  Spectra  of  22.7  kg  TNT  Explosion  ( I )  at  the 
Reference  Range  dim),  (2)  after  Propagating  23  km 
Neglecting  Finite  Amplitude  Effects  (Case  A),  and  (3) 
after  Propagating  23  km  Considering  Finite  Amplitude 
Effects  (Case  D) 


6;s»l 


-*!1 
.l:r  ' 


H.!  ■  i 

>.t.f 

# 

m 

■Sr‘ 


m 


•*U! 

15 

ij 

ij 

‘•5 

•w 


>*■1 


■V* 

I*. 


6. 1 4  Energy  Spectra  of  22.7  kg  TNT  Explosion  after  Propagating 
23  km  Considering  Finite  Amplitude  Effects  ( 1 )  over  the 
First  150  m  (Case  B)  and  (2)  over  the  Entire  23  km  (Case  D) 

6. 1 5  Energy  Spectra  of  22.7  kg  TNT  Explosion  after  Propagating 
23  km  Considering  Finite  Amplitude  Effects  over  ( 1 )  the 
First  1 100  m  (Case  C)  and  (2)  the  Entire  Propagation 

Path  (Case  D) 

A 1  The  Derivatives  of  the  Ray  Path  Position  Vector 
B.l  The  Nonlinearity  of  the  Pressure  Density  Relationship 


I  ’Vi 


I 


p's  P  -  P, 


r 

^0 

r 

R 

R« 


s 


f  s  t  -  t(r) 


T 

T 


u 

W 

X 

z 

Z 

a  ■  p(PoCo®) 


P 

Ps 


Static  pressure 
Acoustic  pressure 
Radial  range 
Reference  range 

Ray  position  vector 
Radius  of  curvature 

Characteristic  distance 
Ray  path  length 
Reference  path  length 
Time 

Retarded  time  variable 
Characteristic  time 

Temperature 
Tangent  to  ray  path 
Fluid  particle  velocity 
Transformed  pressure 
Distance  in  x  direction 
Depth 

Distortion  range  variable 
Convenient  thermodynamic  variable 
Coefficient  of  nonlinearity 
p  at  source  position 


n 


Kronecker  delta 

$ 

Acute  angle  between  ray  path  and 
horizontal 

K 

Thermal  conductivity 

X 

Dllatatlonal  viscosity  coefficient 

Shear  viscosity  coefficient 

Salinity 

P 

Total  density 

Po 

Static  density 

P'»p-Po 

Acoustic  density 

a 

Nondimenslonal  distance 

I 

Viscous  stress  tensor 

♦  ,  t 

Ray  launch  angles 

X 

Entropy 

t(r) 

Elkonal 

FOREWORD 


This  report  is  adapted  from  the  master's  thesis  of  the  same  title 
by  Frederick  D.  Cotaras.  Beginning  in  September  1983  Mr.  Cotaras  was 
enrolled  in  the  Department  of  Electrical  and  Computer  Engineering.  He 
received  his  degree  in  August  1985.  During  this  two-year  period  he  was 
on  educational  leave  from  Defence  Research  Establishment  Atlantic 
(DREA),  Halifax.  N.  S..  Canada,  and  was  supported  Jointly  by  DREA  and  a 
Canadian  Natural  Science  and  Engineering  Research  Council  Postgraduate 
Scholarship. 

The  research  was  carried  out  at  Applied  Research  Laboratories  and 
was  supported  in  part  by  the  Office  of  Naval  Research  (ONR)  under 
Contracts  N00014-75-C-0867  and  N00014-84-K-0574.  The  Scientific 
Officer  for  ONR  was  L.  E.  Hargrove. 

David  T.  Blackstock 
Supervisor 


CHAPTER  1 


INTRODUCTION 

The  subject  of  this  report  is  the  propagation  of  intense  acoustic 
signals  through  an  inhomogeneous  ocean.  The  analysis  is  conducted  both 
analytically  and  numerically.  Our  approach  is  to  use  nonlinear  geometrical 
acoustics.  The  self-refraction  of  acoustic  waves  is  neglected,  and  the 
acoustic  field  is  assumed  to  consist  of  outgoing  waves  only.  The  effects  of 
reflections  from  the  ocean  bottom  and  surface  and  the  effects  of  the 
focusing  of  acoustic  rays  near  caustics  are  not  considered.  In  our 
numerical  algorithm  we  assume  a  stratified  ocean  even  though  the  analytic 
development  is  for  a  fully  inhomogeneous  ocean  in  which  temperature, 
salinity,  density,  and  sound  speed  vary  with  position.  No  losses  are 
considered  in  the  analytic  work.  In  the  numerical  algorithm,  however,  the 
losses  due  to  viscosity  and  relaxation  are  accounted  for  in  an  ad  hoc 
fashion.  Our  algorithm  is  designed  to  propagate  signals,  which  may  contain 
weak  discontinuities,  along  a  ray  path  using  a  stepwise  time  domain 
technique  and  to  apply  absorption  in  the  frequency  domain. 

When  considering  how  to  approach  a  problem  involving  an  intense 
sound  wave,  one  usually  examines  the  corresponding  small-signal  problem. 


Popular  techniques  for  dealing  with  the  propagation  of  small-signal  sounds 
Include  the  following:  normal  modes,  parabolic  approximation  to  the  wave 
equation,  and  geometrical  acoustics.  Since  we  are  Interested  In  handling 
discontinuities,  an  approach  which  lends  itself  to  a  time  domain  implemen¬ 
tation  is  preferred.  Thus  geometrical  acoustics  and  the  parabolic  approxi¬ 
mation  are  both  good  avenues  for  dealing  wih  intense  sounds.  As  noted 
above,  we  have  chosen  to  use  geometrical  acoustics.  The  alternative 
approach,  the  parabolic  approximation,  is  being  explored  by  McDonald  and 
Kuperman(1984). 

Our  implementation  of  nonlinear  geometrical  acoustics  is 
somewhat  restricted.  Since  self-refraction  is  neglected,  the  ray  paths  are 
assumed  to  be  determined  solely  by  the  inhomogeneous  medium  and  not 
affected  by  the  amplitude  of  the  signal  (Whitham  1956).  It  turns  out  that 
this  assumption  places  a  limit  of  282  dB// 1  pPa  on  the  peak  amplitude  of 
the  signal.  Reflections  and  focusing  are  neglected  because  neither  are  fully 
understood  for  the  case  of  intense  sound  waves.  Reflections  and  focusing 
are  avoided  in  the  numerical  work  by  a  careful  selection  of  the  acoustic  ray 
paths.  Because  the  acoustic  rays  are  not  permitted  to  either  interact  with 
ocean  surface  and  bottom  or  to  pass  through  caustics,  the  propagation 
range  is  limited  to  a  maximum  of  about  70  km  in  the  deep  ocean. 

The  primary  objective  of  this  work  is  to  investigate  the 
importance  of  nonlinear  effects  at  long  ranges  from  an  intense  acoustic 
source  such  as  an  explosive.  Because  the  amplitude  of  most  underwater 
sounds  is  small,  the  most  commonly  used  theory  in  underwater  acoustics  is 
small-signal  theory.  When  dealing  with  underwater  explosions. 


Investigators  commonly  assume  that,  after  propagating  a  certain  distance, 
the  Intense  sound  Is  sufficiently  diminished  that  small-signal  theory 
applies.  One  of  the  goals  of  this  thesis  Is  to  attempt  to  quantify  this 
assumption. 

A.  Review  of  Finite  Amplitude  Effects 

An  Intense,  or  finite  arr^Htucte,  acoustic  signal  Is  distinguished 
from  the  better  known  small,  or  infinitesimal  signal  by  the  amplitude  of 
the  field  variables  such  as  the  pressure  or  particle  velocity.  To  aid  In  the 
distinction,  we  Introduce  the  acoustic  Mach  number  M,  a  nondimenslonal 
number  equal  to  the  ratio  of  the  maximum  fluid  particle  velocity  to  the 
small-signal  sound  speed.  In  the  case  of  small-signals,  M  Is  approximately 
zero,  infinitesimally  small.  However,  for  the  finite  amplitude  signals 
considered  In  this  report,  M  may  be  as  large  as  0.06,  a  small,  but  finite 
value. 

It  Is  well  known  that  the  propagation  of  a  finite  amplitude  wave 
cannot  be  accurately  modeled  by  a  small-signal  acoustical  theory.  The 
reason  small-signal  theory  fails  Is  that  It  does  not  correctly  predict  the 
propagation  speed  of  the  wave.  The  physical  mechanisms  that  cause  the 
propagation  speed  to  vary  from  the  small-signal  value  are  ( I )  the 
nonlinearity  of  the  pressure-density  relationship  of  the  medium  and  (2) 
convection.  Convection  occurs  when  the  fluid  particles  themselves  are  set 
into  motion  by  the  passing  acoustic  wave  and  contribute  their  velocity  to 
the  total  wave  speed.  The  two  phenomena,  the  nonlinearity  of  the 
pressure-density  relationship  and  convection,  are  embodied  in  the 


coefficient  of  nonlinearity  p.  These  phenomena  are  present  even  In  the 
cases  of  small  signals;  they  are  just  not  perceivable.  As  the  amplitude  of 
the  signal  Increases,  so  does  the  significance  of  their  effects.  To 
incorporate  these  effects  Into  acoustical  theory,  one  must  consider  the 
nonlinear  terms  In  the  hydrodynamics  equations.  Therefore  the  theory  that 
governs  a  finite  amplitude  wave  must  be  a  nonlinear  theory.  Hence  the 
words  finite  amplitude  and  nonlinear  are  used  Interchangeably  In  this 
report.  Similarly  the  words  small-signal  and  linear  are  used 
Interchangeably  since  small-signal  theories  are  derived  neglecting  all 
terms  except  linear  terms. 

The  problem  of  propagation  of  plane  finite  amplitude  waves 
through  a  homogeneous  medium  has  a  solution  that  Is  well  established;  see, 
for  example,  Blackstock  (1972).  The  nonlinear  acoustical  theory  described 
by  Blackstock  Is  referred  to  as  weak-shock  theory.  A  shock  is  a 
discontinuity  In  the  field  variable.  If  the  shock  is  small  enough  to  be  dealt 
with  by  a  nonlinear  theory  that  is  quadratic  In  the  field  variable,  the  shock 
Is  referred  to  as  weak  (see,  for  example,  Whitham  1974,  p.  37).  The  highest 
amplitude  dealt  with  In  this  report  corresponds  to  M  =  0.06  or, 
equivalently,  a  peak  pressure  level  of  282.6  dB  // 1  pPa.  According  to 
Pestorlus  and  Williams  (1974),  this  level  Is  within  the  amplitude  limits  of 
weak-shock  theory.  Pestorlus  (1973)  developed  a  computer  based  version 
of  the  weak-shock  solution  to  plane  finite  amplitude  acoustic  propagation. 
He  used  the  computer  program  to  solve  problems  Involving  both  finite 
amplitude  noise  and  periodic  waves.  As  previously  noted,  a  modified  form 
of  Pestorlus's  algorithm  Is  used  in  this  report. 
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The  extension  of  finite  amplitude  plane  wave  theory  to  nonplanar 
waves  Is  of  notable  Importance  to  this  work  (see,  for  example,  Blackstock 
1964).  Blackstock  extended  the  plane  wave  solution  to  the  problem  of 
spherical  and  cylindrical  finite  amplitude  acoustic  waves.  His  approach 
was  to  use  two  transformations  to  reduce  the  spherical  (cylindrical)  wave 
equation  to  the  form  of  the  finite  amplitude  plane  wave  equation.  The  two 
transformations,  one  for  the  field  varl2tf)le  and  the  other  for  the  range, 
permit  the  use  of  the  solution  of  the  plane  wave  problem.  The 
transformations  required  to  deal  with  plane  waves  propagating  vertically 
In  a  stratified  medium  were  developed  by  Carlton  and  Blackstock  (1974). 

The  Idea  of  reducing  a  complex  problem  to  a  simpler  one  with  a  known 
solution  Is  used  in  this  thesis. 

B.  Review  of  Linear  Geometrical  Acoustics 

Two  separate  geometrical  acoustic  theories,  one  for  infinitesimal 
signals  and  one  for  finite  amplitude  signals,  exist.  The  extension  of  linear 
geometrical  acoustics  to  Incorporate  finite  amplitude  signals  Is  a 
comparatively  recent  development.  The  two  theories  predict  different 
wave  shapes,  but  the  same  ray  paths.  The  older  of  the  two,  linear 
geometrical  acoustics.  Is  described  in  the  literature  survey  that  follows. 

Several  methods  for  developing  linear  geometrical  acoustics 
exist.  One  approach  originates  in  geometrical,  or  short-wave,  optics  as 
developed  by  Hamilton  ( 1 832;  see  Conway  and  Synge  1931).  The  central 
equation  In  the  geometrical  technique  for  both  acoustics  and  optics  Is  the 
so  called  elkonal  equation.  The  elkonal  equation  defines  the  ray  paths  (see. 


for  example,  Bom  and  Wolf  1980,  p.  1 12).  The  word  eikonal  comes  fom  the 
Greek  word  eikon,  for  Image,  and  was  Introduced  Into  physics  by  Bruns 
(1895)  (see,  for  example,  Sommerfeld  1949,  p.  207).  Bruns  developed  a 
function  similar  to  the  one  developed  by  Hamilton  and  called  It  the  eikonal. 
The  eikonal  Is  a  function  of  the  spatial  coordinate  system  and  defines  "a 
system  of  surfaces  the  orthogonal  trajectories  of  which  are  rays" 
(Sommerfeld  1 949,  p.  337).  Sommerfeld  and  Runge  (1911)  were  the  first  to 
derive  the  eikonal  equation  from  the  scalar  wave  equation.  Their  derivation 
clearly  shows  that  ray  theory  Is  exact  only  In  the  case  of  Infinite 
frequency.  Ray  theory  Is,  however,  a  valid  approximation  to  surprisingly 
low  frequencies.  It  turns  out  that  the  approximation  depends  on  the  degree 
of  Inhomogeneity  In  the  medium. 

Sommerfeld  (1949,  p.  210)  provides  some  Interesting  comments 
about  the  eikonal  and  linear  geometrical  techniques  In  general.  He  states 
that  the  eikonal  equation  Is  a  specialization  of  Hamilton's  equation  of 
dynamics.  He  notes  that  Hamilton  worked  on  optical  problems,  then  went  on 
to  apply  his  knowledge  to  dynamics.  Hamilton's  work  In  dynamics  Is 
embodied  In  the  Hamllton-Jacobi  equation,  and  analogies  between  the 
eikonal  and  the  action  of  the  material  particle  are  often  made.  Goldstein 
( 1 950,  p.  3 1 2)  states  that  the  Hamilton- Jacobi  equation  "tells  us  that 
classical  mechanics  corresponds  to  the  geometrical  optics  limit  of  a  wave 
motion."  Consequently  geometrical  acoustics  Is  sometimes  derived  via  the 
Hamilton- Jacobi  equation  (see,  for  example,  Landau  and  Lifshitz  1959, 
p.  257).  Sommerfeld  observes  that  the  WKB  (Wentzel-  Kramers-  Brillouin) 


approximation  of  the  wave  equation  corresponds  to  a  transition  from  wave 
optics  (wave  acoustics)  to  geometrical  optics  (geometrical  acoustics). 

Another  way  to  develop  linear  ray  theory  is  to  use  Huygens' 
principle  and  Snell's  law.  The  first  application  (within  acoustics)  of  this 
approach  was  to  analyze  the  problem  of  propagation  in  a  moving, 
inhomogeneous  atmosphere  (see,  for  example.  Stuff  1979).  The  following  is 
a  brief  summary  of  some  derivations  and  applications  of  linear  geometrical 
acoustics  developed  in  this  fashion. 

Rayleigh  (1896,S  289)  examined  the  problem  of  acoustic 
propagation  in  a  windy  atmosphere.  However  he  did  not  acknowledge  the 
difference  between  the  direction  of  the  acoustic  ray  path  and  that  of  the 
normal  to  the  wave  (eikonal)  surface.  This  flaw  has  been  commented  on  by 
several  authors  (Barton  1901;  Kornhauser  1953;  Lighthill  1965;  Thompson 
1972;  Ugincious  1972  ).  Several  other  researchers  made  use  of  linear  ray 
theory  to  study  propagation  of  airborne  sound.  An  early  application  is  given 
by  Fujiwhara  ( 1 9 1 2;  1 9 1 6).  He  derived  three-dimensional  acoustic  ray 
theory  and  applied  it  to  the  problem  of  sounds  produced  by  the  volcano 
Asama  in  central  Japan.  Another  interesting  early  work  is  that  of  Milne 
( 1921 ).  His  work  is  an  outgrowth  of  research  done  during  the  First  World 
War  on  the  acoustic  detection  of  aircraft.  Rothwell  (1947)  applied  acoustic 
ray  theory  to  the  problem  of  meteorological  investigations  by  acoustical 
techniques.  It  is  interesting  to  note  that  Rothwell's  data  is  from 
experiments  performed  in  1930  and  1931  in  which  16  lb  practice  shells 
were  fired  from  an  anti-aircraft  gun.  The  data  were  originally  applied  to 
the  problem  of  acoustical  detection  of  aircraft. 


Because  of  the  need  to  detect  submarines,  research  on  linear 
geometrical  acoustics  underwent  a  resurgence  during  the  Second  World 
War.  In  the  years  immediately  following  the  war,  much  of  this  research 
was  published.  In  a  major  work  Blokhintzev  (1946a;  1946b)  derived  the 
elkonal  equation  for  a  moving,  inhomogeneous  medium.  Earlier  researchers 
had  assumed  that  the  acoustic  energy  is  constant  along  the  ray  path. 
Blokhinstzev  found  that  a  certain  combination  of  the  acoustic  pressure,  ray 
tube  area,  sound  speed,  particle  velocity,  and  density  (not  necessarily  the 
energy)  stays  constant  along  the  ray  path.  From  his  work  the  term 
Blokhintzev  invariant  arose.  Frank,  Bergmann,  and  Yaspen  (1969)  discussed 
the  derivation  of  ray  theory  and  gave  a  good  account  of  the  frequency  limits 
of  the  theory.  Bergmann  ( 1 946)  discussed  the  Importance  of  the  density 
variation  and  the  gravity  terms  In  the  derivation  of  ray  theory. 

Over  the  years  since  the  end  of  the  Second  World  War,  several 
other  derivations  of  linear  ray  theory  have  been  given.  Using  tensor 
calculus  and  a  general  curvilinear  coordinate  system,  Haselgrove  (1954) 
derived  ray  theory  via  reciprocal  surfaces.  He  applied  his  theory  to  radio 
wave  propagation,  specifically  to  the  calculation  of  ionospheric  ray  paths,  a 
requirement  for  determining  the  maximum  usable  frequencies  for 
short-wave  radio  transmission.  Haselgrove  appears  to  have  been  the  first 
to  use  an  electronic  computer  (EDSAC,  Cambridge  University  Mathematical 
Library)  to  calculate  ray  paths. 

More  recent  derivations  of  ray  theory  Include  those  of  Eby  and 
Mal'tsev.  Eby  (1967)  derived  three-dimensional  ray  tracing  using  a  Frenet 
formulation.  In  another  paper  (Eby  1970)  he  examined  the  ray  paths  as 


temporal  geodesics.  Martsev  (1983)  derived  equations  In  barycentric 
coordinates,  a  technique  that  simplifies  the  calculation  of  the  ray  paths. 


C.  Geometrical  Acoustics  for  Finite  Amplitude  Signals 

Until  now  the  discussion  has  centered  on  the  problem  of  deriving 
ray  theory  for  small-signal  waves.  We  now  address  the  problem  of 
geometrical  acoustics  for  finite  amplitude  signals.  Heller  (1953)  derived 
the  equation  for  the  ray  path  of  an  acoustic  discontinuity  propagating  In  a 
moving,  Inhomogeneous  medium.  He  referred  to  this  equation  as  a 
"generalized  elkonal  equation"  and  showed  that  weak  shocks  move  along  the 
same  ray  paths  as  small-signal  acoustic  waves.  Komhauser  (1953)  showed 
that  Heller's  generalized  elkonal  equation  could  be  obtained  by  a  simple 
extension  of  the  small-signal  elkonal  equation.  Keller  (1954)  was  the  first 
to  discuss  the  variation  of  the  shock  strength  along  the  ray  path.  Keller 
solved  the  problem  by  examining  the  Jump  of  the  field  variables  across  the 
surface  of  discontinuity  (see,  for  example,  Jeffrey  and  Tanluti  1964).  This 
approach  Is  not  common  In  acoustics  today.  Whitham  (1956),  In  a  major 
work  on  the  propagation  of  weak  shocks,  described  extensions  of  linear 
geometrical  acoustics  to  encompass  weak  shocks.  His  development  Is  very 
closely  related  to  the  nonlinear  geometrical  acoustical  theory  developed  In 
this  report.  Whitham  assumed,  as  we  do,  that  the  self-refraction  of  the 
acoustic  ray  paths  Is  Insignificant.  Thus  he  used  the  ray  paths  from  linear 
ray  theory  and  went  on  to  develop  a  set  of  Improved  relations  to  predict  the 
shock  strength.  In  a  discussion  on  the  Inadequacies  of  linear  geometrical 
acoustics,  Whitham  said,  "It  should  be  stressed  that  this  Inaccuracy  Is  a 
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failure  of  the  linear  theory  of  sound  and  Is  not  Introduced  by  the 
approximations  of  geometrical  acoustics.” 

Also  using  a  geometrical  technique,  GiMln  (1958)  addressed  the 
problem  of  propagation  of  acoustic  discontinuities  In  an  Inhomogeneous 
medium.  The  solution  was  presented  In  a  way  more  familiar  to  researchers 
In  acoustics,  that  Is,  by  starting  with  the  hydrodynamical  equations  and 
making  simplifying  assumptions  based  on  the  degree  of  smallness  of  the 
terms  In  the  equations.  Using  similar  notation,  Ostrovsky  ( 1 963)  examined 
linear  geometrical  acoustics  In  nonstationary  media.  He  then  extended  the 
theory  to  Include  finite  amplitude  waves  and  applied  the  theory  to  the 
problem  of  a  sinusoidal  finite  amplitude  wave  In  a  steadily  moving  gas. 
Ostrovsky  Introduced  a  parameter  ”g"  which  Indicates  the  Influence  of 
medium  Inhomogeneity  upon  the  shock  formation  distance. 

In  the  late  I960's,  there  was  a  large  Interest  In  the  propagation  of 
finite  amplitude  acoustic  waves  In  an  Inhomogeneous  moving  atmosphere. 
The  stimulation  was  the  problem  of  sonic  booms  from  supersonic 
transports  (SST).  Hayes,  Haefell,  and  Kulsrud  ( 1 969)  developed  a  computer 
program  that  uses  small-signal  acoustic  ray  theory,  and  then  modifies  the 
result  to  account  for  finite  amplitude  effects.  Seebass  (1969)  presented 
the  derivation  of  the  equation  for  the  acoustic  pressure  behind  a 
propagating  weak  shock.  The  equation  was  obtained  by  correcting  small- 
signal  ray  theory  results  for  nonlinear  effects.  Seebass's  result  shows 
how  the  acoustic  pressure  depends  on  the  acoustic  ray  tube  area  In  the  case 
of  finite  amplitude  waves. 


In  the  mid-1970’s,  interest  in  finite  amplitude  propagation  in 
stationary  inhomogeneous  media  reappeared.  Carlton  and  Blackstock  (1974) 
analyzed  the  problem  of  vertical  propagation  of  finite  amplitude  plane 
waves  in  a  horizontally  stratified  ocean  The  transform  technique  that 
Carlton  and  Blackstock  used  is  similar  to  that  used  by  Blackstock  (1964)  for 
nonplanar  waves  in  a  homogeneous  medium.  In  a  comment  in  the  abstract, 
Carlton  and  Blackstock  (1974)  state  that  their  results  were  intended  to  be 
incorporated  in  a  ray  theory  for  finite  amplitude  acoustic  waves  in  the 
inhomogeneous  ocean. 

Before  Carlton  and  Blackstock  had  an  opportunity  to  perform  this 
extension,  however,  the  results  of  a  parallel  development  were  disclosed 
Ostrovsky,  Pelinovsky,  and  Fridman  (1975;  1976)  described  the  solution  of 
the  problem  of  propagation  of  finite  amplitude  waves  in  a  stationary, 
inhomogeneous  medium.  Their  work  is  central  to  the  discussion  of 
nonlinear  geometrical  acoustics  (N6A)  in  this  report. 

Nonlinear  geometrical  techiques  have  been  used  in  many 
investigations.  Ostrovsky  (1976)  discussed  applications  of  nonlinear 
geometrical  techniques  to  problems  such  as  acoustic  propagation  in 
nonstationary  media  and  the  heating  of  the  sun  chromosphere.  Pelinovsky, 
Petukhov,  and  Fridman  (1979)  extended  NGA  to  include  the  effects  of  ocean 
salinity.  NGA  was  also  expanded  to  encompass  the  effects  of  small 
dissipation  and  dispersion  (Pelinovsky  and  Fridman  1983).  Warshaw  (1980) 
extended  the  results  of  Blokhintzev  ( 1 946a,  1 946b)  to  account  for  the 
cumulative  effects  of  the  dissipative  and  second-order  convective  terms. 
More  recently  the  work  of  Pelinovsky,  Petukhov,  and  Fridman  (1979)  was 


examined  by  Morfey  ( t984a).  Morfey  Introduced  a  parameter  'G"  to  describe 
the  effect  of  the  Inhomogeneity  on  the  shock  formation  distance.  He 
calculated  the  parameter  for  a  variety  of  oceanic  environments  and 
concluded  that  the  effect  of  Inhomogeneity  upon  finite  amplitude 
propagation  Is  small.  Morfey's  work  Is  used  In  this  report. 

0.  Scope  of  the  Investigation 

This  report  is  divided  as  follows:  In  Chapter  2  we  establish  the 
criteria  for  simplifying  the  equations  of  hydrodynamics.  The  equations  are 
then  simplified  for  the  cases  of  small  signals  In  homogeneous  fluids,  small 
signals  in  inhomogeneous  fluids,  and  finite  amplitude  signals  in 
inhomogeneous  fluids.  In  Chapter  5  we  discuss  linear  geometrical 
acoustics  and  derive  the  eikonal  and  transport  equations.  At  the  same  time 
we  introduce  some  of  the  techiques  required  for  the  case  in  which  the 
signals  are  of  finite  amplitude.  Nonlinear  geometrical  acoustics  is 
presented  in  Chapter  4.  The  notation  used  varies  slightly  from  that  of 
Pelinovsky  et  al.  (1979)  and  adheres  more  closely  to  that  of  Blackstock 
(1964).  In  Chapter  5  we  discuss  the  numerical  implementation  of  the 
findings  of  Chapter  4.  The  method  of  handling  shocks  in  the  waveform  and 
the  transformation  to  the  frequency  domain  to  correct  for  absorption  is 
due  to  Pestortus  ( 1 973).  The  computer  based  ray  model  used  is  due  to 
Foreman  ( 1 983)  Some  of  the  results  of  the  testing  of  the  numerical 
algorithm  are  presented.  In  Chapter  6  the  numerical  algorithm  is  used  to 
study  the  effects  of  viscosity,  relaxation,  and  dispersion  on  the  propagation 
of  waves  through  the  ocean.  It  is  also  used  to  provide  examples  of  the 
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effects  of  ocean  Inhomogeneity  on  finite  amplitude  propagation.  The 
waveforms  considered  are  a  weak  shock  with  an  exponentially  decaying  tail 
and  Morfey's  (1985)  modified  form  of  the  Wakeley  explosion  waveform 
(Wakeley  1977). 


CHAPTER  2 
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RANKING  OF  TERMS 

A  Introduction 

In  this  chapter  a  method  for  simplifying  the  equations  of 
hydrodynamics  Is  discussed.  An  equation  of  state  for  seawater  and  the 
hydrodynamics  equations,  including  loss  terms,  are  presented.  As 
mentioned  In  Chapter  I,  it  is  our  intention  to  neglect  all  losses  during  the 
analytic  development,  and  then  account  for  viscosity  and  relaxation  by  a 
special  procedure  In  the  numerical  propagation  routine.  The  viscosity  and 
heat  conduction  terms  have,  however,  been  included  in  the  hydrodynamics 
equations  so  that  the  ranking  system  can  be  fully  demonstrated.  Basic 
assumptions  about  ( I )  the  amplitude  of  the  signal,  (2)  the  magnitude  of  the 
loss  coefficients,  and  (3)  the  type  and  degree  of  inhomogeneity  of  the 
medium  are  made.  The  terms  in  the  equations  are  then  ranked  according  to 
their  relative  importance.  Use  of  the  ranking  system  enables  simplified 
forms  of  the  hydrodynamics  equations  to  be  readily  obtained  (Lighthill 
1 956;  Carlton  and  Blackstock  1 974). 

For  an  Inhomogeneous,  thermally  conducting,  viscous  fluid,  an 
equation  of  state  and  the  hydrodynamics  equations  are  as  follows  (see,  for 
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example,  Panton  1984;  the  momentum  equation  is  from  Hunt  1955); 

State  P  =  Pa.X.p)  .  (A.1) 

Continuity  ap/dt  +  V*pu  =  0  ,  (A.2) 

Momentum  pau/at  +  p(u-V)u  *  -VP  ♦  pg  ♦  (X+2p)V(V-u) 

-  pVx(Vxu)  +  (V*u)VX 
♦  2(V|l  V)u  +  Viix(Vxu)  ,  (A.3) 

Energy  pae/at  ♦  p(u*V)e  -  -P(V*u)  ♦  t:Vu  +  V-(KVT)  ,  (A.4) 

where  P  is  the  pressure,  I  is  the  salinity,  X  is  the  entropy,  p  is  the  density, 
t  is  the  time,  u  is  the  particle  velocity,  g  is  gravity,  X  and  p  are, 
respectively,  the  dilatational  and  shear  viscosity  coefficients,  e  is  the 
Internal  energy,  i  is  the  viscous  stress  tensor,  K  is  the  coefficient  of 
thermal  conductivity,  and  T  is  the  temperature.  Since  the  emphasis  of  this 
work  is  on  acoustic  propagation  in  the  ocean,  the  salinity  of  the  ocean  Is 
included  In  the  equation  of  state.  Note  also  that  our  coordinate  system  has 
gravity  acting  in  the  positive  z  direction;  that  is,  z  is  positive  downward. 

We  expand  the  continuity  and  momentum  equations.  First  we 
express  the  density  as  the  sum  of  a  static  value  Pq  and  a  small  fluctuation 
p',  and  the  pressure  as  the  sum  of  a  static  value  Pq  and  a  small  fluctuation 
P', 


psPo^P’ 


(A5) 


P  ■  pQ  ♦  P' 


(A6) 


In  a  static  fluid  the  pressure  and  density  fluctuations  are  zero.  In  this 
case,  the  continuity  equation,  Eq.  (A2),  is  identically  satisfied  while  the 
momentum  equation,  Eq.  (A3),  reduces  to 

VPo  '  PoO  •  (A.7) 

If  Eqs.  (A5)  and  (A6)  are  substituted  into  the  continuity  and  momentum 
equations,  and  Eq.  (A7)  is  used  to  reduce  terms  in  the  latter,  the  following 
equations  are  obtained: 

Continuity  dp'/dt  ♦  p^V-u  *  -u-Vp^  -  u*Vp'  -  p' V*u  ,  (A8) 

Momentum  p^au/at  ♦  VP'  -  (p'/Po)VPo  -  p’du/at  -  Po(u*V)u  -  p'(u*V)u 

♦  (X^2|i)V(V-u)  -  pVx(Vxu) 

♦  (V*u)VX  ♦  2(Vp'V)u  ♦  Vpx(Vxu) 

♦  (V-u)V(X  ♦  2p)  .  (A9) 

To  simplify  Eqs.  (A8)  and  (A9),  we  must  have  some  knowledge  of  the 
relative  importance  of  each  of  the  terms. 

We  now  state  our  basic  assumptions;  The  magnitude  of  the 
particle  velocity,  M,  is  assumed  to  be  small  with  respect  to  the  sound 
speed  Cq,  but  not  necessarily  infinitesimal.  It  turns  out  that  this 
assumption  implies  that  the  pressure  fluctuation  P'  is  small  with  respect 
to  PqOo^  3nd  that  the  fluctuation  p’  is  small  with  respect  to  Pq.  The  above 
assumption  is  expressed  mathematically  as  follows: 


(A.11) 
(A.  1 2) 


Ip'UPo  . 

IP‘l«PoCo^  • 

As  for  the  inhomogeneity  of  the  medium,  it  is  assumed  that  environmental 
properties  of  the  seawater  such  as  the  density,  temperature,  salinity,  and 
sound  speed  vary  with  positioa  As  can  be  seen  from  Eq.  (A.7),  the  static 
pressure  is  assumed  to  vary  only  with  depth.  We  assume  that  the  ocean  is 
only  mildly  inhomogeneous;  that  is,  the  environmental  properties  vary 
slowly  on  a  wavelength  scale.  Consequently  the  derivative  of  an 
environmental  property,  for  example  Vp^,  is  small;  that  is, 

MVPoUPo  .  (A  1 3) 

where  X  is  the  wavelength.  Losses  due  to  heat  conduction  and  the  diffusion 
of  dissolved  salts  are  assumed  to  be  zero;  however,  losses  due  to  viscosity 
and  relaxation  are  assumed  to  be  small,  but  not  zero. 

Before  stating  the  ranking  system,  we  define  some  nomenclature. 
A  linear  term  with  a  coefficient  containing  a  first  derivative  of  an 
environmental  parameter  is  referred  to  as  an  inhomogeneity  term; 
examples  are  u- Vpg  or  VPqCpVpq).  Linear  terms  with  coefficients 
containing  a  constant  loss  coefficient,  such  as  (X  ♦  2p)(V’U),  are  referred 
to  as  dissipation  terms.  Similarly  terms  containing  a  quadratic  nonlinear 
term  with  a  coefficient  that  contains  neither  a  derivative  of  an 
environmental  property  of  the  fluid  nor  a  loss  coefficient  are  referred  to 


as  nonlinearity  terms;  examples  are  Pq(u-V)u  and  p'(V’U).  Nonlinearity 
terms  are  associated  with  finite-amplitude  behavior. 

Use  of  the  assumptions  stated  above  enables  the  various  terms  In 
the  hydrodynamics  equations  to  be  ranked  by  degree  of  smallness  as 
first-order  terms,  second-order  terms,  and  higher-order  terms.  A 
first-order  term  Is  defined  as  a  linear  term  with  a  coefficient  that 
Involves  neither  derivatives  of  environmental  properties  of  the  fluid  nor  a 
loss  coefficient.  Examples  of  first-order  terms  are  Po(V*u)  and  pQdu/dt. 
Second-order  terms  represent  only  the  most  Important  effects  of 
nonlinearity.  Inhomogeneity,  and  losses.  Accordingly  Inhomogeneity  terms, 
nonlinearity  terms,  and  dissipation  terms  are  classified  as  second-order 
terms.  Since  the  effects  of  nonlinearity.  Inhomogeneity,  and  losses  were 
assumed  small,  any  term  representing  the  Interaction  of  any  two  effects 
would  be  expected  to  be  negligible.  Such  terms  are  encompassed  In  the 
third  major  category,  higher-order  terms.  Examples  of  higher-order  terms 
are  cubic  terms  such  as  p'(U'V)u  and  nonlinear  Inhomogeneity  terms  such 
as  VPq(u*V)u.  Linear  terms  with  coefficients  that  contain  higher-order 
derivatives  of  environmental  properties,  or  that  are  nonlinear  In  the  first 
derivative  of  an  environmental  property  are  also  higher-order  terms; 
examples  are  (VPq*VPq)p'  and  p'V^Pq. 

Use  of  the  ranking  system  makes  It  relatively  easy  to  deal  with 
special  cases  such  as  small-signal  waves  In  a  lossless  Inhomogeneous  fluid 
and  finite  amplitude  (but  not  strong)  waves  in  dissipative  homogeneous 
fluid.  The  hydrodynamics  equations  can  be  simplified  and.  In  most  cases, 
combined  to  form  the  wave  equation  corresponding  to  the  situation  at  hand 


The  analysis  from  this  point  on  deals  solely  with  lossless  fluids. 
The  effects  of  the  dissipation  terms  are  considered  In  the  numerical 
propagation  routine  discussed  In  Chapter  5.  In  this  section  the 
hydrodynamics  equations  are  simplified  to  find  ( I )  linear  equations  for  a 
homogeneous  fluid,  (2)  linear  equations  for  an  Inhomogeneous  fluid,  and 
(3)  nonlinear  equations  for  an  Inhomogeneous  fluid.  Since  all  three  sets  of 
equations  are  for  lossless  fluids,  the  loss  terms  can  be  dropped  from  the 
expanded  momentum  equation,  Eq.  (A.9),  and  from  the  energy  equation, 

Eq.  (A.4).  The  state  and  hydrodynamics  equations  are  then  as  summarized 


below; 

State  P»P(t,X,p)  ,  (A.I) 

Continuity  dp'/dt  ♦  Pq(V*u)  »  -u-Vp^  -  u-Vp’  -  p‘(V-u)  ,  (A.8) 

Momentum  p^du/dt  ♦  VP'  -  (p'/Pq)VPo  -  p'du/dt  -  Pq(u-V)u  ,  (B.l) 

Energy  DX/Dt  -  0  ,  (B.2) 

where  D(-)/Dt  Is  the  material  derivative  and  Is  defined  as  follows; 

D(  )/Dt  -  d(  )/dt  ♦  (U-V)(  )  .  (B.3) 


A  physical  Interpretation  of  the  energy  equation  for  an  Inhomogeneous 
fluid,  Eq.  (B.2),  is  that  the  entropy  of  any  given  fluid  particle  remains 
constant  (see,  for  example.  Pierce  1961,  p.  1 2).  In  the  case  of  a 
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homogeneous  fluid,  the  energy  equation  would  be  a  little  simpler.  It  would 
be  Interpreted  as  meaning  that  the  entropy  Is  the  same  for  all  fluid 
particles. 

Since  there  are  five  variables  and  only  four  equations,  we  need 
another  equation.  We  can  get  It  from  our  earlier  assumption  that  the  losses 
due  to  the  diffusion  of  dissolved  salts  are  zero.  This  assumption  Is 
equivalent  to  saying  that  the  change  In  salinity  of  any  one  fluid  particle  Is 
zero  (Landau  and  Lifshitz  1959,  sec.  57), 

Dim  =  0  .  (B.4) 

1.  Linear  Hydrodynamics  Equations  for  Lossless  Homogeneous 

Retention  of  only  first-order  terms  simplifies  the  hydrodynamics 
equations  Into  a  set  of  linear  equations  for  homogeneous  fluids.  Since  all 
the  terms  on  the  right-hand  side  of  Eqs.  (A.8)  and  (B.  1 )  are  second-order 
terms,  they  must  be  dropped.  This  leads  to  the  following: 


Linear  Continuity 

apVdt  ♦  PqV-u  =  0  , 

(B.5) 

Linear  Momentum 

p^du/dt  ♦  VP’  *  0  . 

(B.6) 

Turning  our  attention  to  the  equation  of  state,  Eq.  (A.  1 ),  we  see  that  by  using 
a  Taylor  series  expansion  and  retaining  terms  up  to  first  order,  the 
following  equation  can  be  obtained: 


The  small-signal  sound  speed  Is  defined  as  follows; 


W 


Ui)*Po 


(B.8) 


Substituting  Eq.  (B.8)  Into  the  equation  of  state,  Eq.  (B.7),  and  using  the 
expressions  for  the  density  and  pressure  given  In  Eqs.  (A.5)  and  (A.6),  we 
arrive  at  the  following: 

P'  =  Co2p'  .  (B.9) 

Equations  (B.5),  (B.6),  and  (B.9)  are  the  linearized  continuity,  momentum, 
and  state  equations  for  a  lossless  homogeneous  ocean. 

2.  Linear  Hydrodynamics  Equations  for  Lossless  Inhomogeneous 

Fluids 

Retention  of  first-order  terms  and  Inhomogeneity  terms 
simplifies  the  hydrodynamics  equations  Into  a  set  of  linear  equations  for 
Inhomogeneous  fluids.  The  continuity  equation  Is  obtained  from  Eq.  (A.8)  by 
dropping  the  nonlinear  terms,  the  last  two  terms  on  the  right-hand  side. 
Similarly,  the  momentum  equation  is  obtained  from  Eq.  (B.  I ): 

Continuity  dpVdt  ♦  p^V-u  ♦  u-Vp^  =  0  ,  (B.  10) 

p^du/dt  ♦  VP*  -  (p7pq)VPq  *  0 


Momentum 


(B.1I) 
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The  equation  of  state  Is  obtained  by  taking  the  material  derivative  of 
Eq.(A.l), 


DP 

Dt 


Dp 

Dt 


ff)  • 


DX 

Dt 


(B.12) 


If  Eqs.  (B.2)  and  (B.4)  are  substituted  into  Eq.  (B.  12),  the  second  and  third 
terms  of  Eq.  (B.  1 2)  drop  out.  Expansion  using  the  definition  of  the 
smali-signal  sound  speed.  Eq.  (B.8).  yields 


State  dP'/dt  ♦  u-VpQ  =  0^2  (dpVdt  ♦  u* Vp^)  .  (B.  1 3) 

Thus  the  linear  equations  of  continuity,  momentum,  and  state  for  a  lossless, 
inhomogeneous  fluid  are  Eqs.  (B.  1 0),  (B.  1 1 ),  and  (B.  1 3). 

3.  Nonlinear  Hudrodunamics  Equations  for  Lossless  Inhomogeneous 
Fluids 

Retention  of  first-order  terms  as  well  as  inhomogeneity  and 
nonlinearity  terms  simplifies  the  hydrodynamics  equations  into  a  set  of 
nonlinear  equations  for  inhomogeneous  fluids.  The  nonlinear  equations  of 
continuity  and  momentum  for  a  lossless  inhomogeneous  fluid  are  Eqs.  (A.8) 
and  (B.  I );  no  terms  need  to  be  dropped.  The  equation  of  state  is  more 
complicated  than  it  was  in  the  previous  two  cases.  Retaining  terms  up  to 
second-order  in  the  Taylor  series  expansion  of  Eq.  (A.1 ),  we  obtain 
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P-P, 


^  ^''X.54>“Po  ''X.^4>“Po 


^  a  il  • 


♦(X-Xq) 


P.X.|*5o 


W  e,  , 

^  P.5.X  -  X, 


(B.14) 


where  It  has  been  assumed  that  variations  In  the  salinity  and  entropy  are  of 
second-order  in  smallness.  Using  the  definition  of  the  small-signal  sound 
speed,  Eq.  (B.8),  and  the  expressions  for  the  density  and  pressure  given  In 
Eqs.  (A.5)  and  (A.6),  we  see  that  Eq.  (B.1 4)  becomes 

Jx.t4>-Po  ^  ^-'p.X^-U 


♦  (X-Xq) 


(B.15) 


The  nonlinear  continuity,  momentum,  and  state  equations  for  a  lossless 
Inhomogeneous  fluid  are  Eqs.  (A.8),  (B.  1 ),  and  (B.  1 5). 


C.  Substitution  Into  Second-Order  Terms  Using  First-Order  Relations 

Once  a  simplified  form  of  the  hydrodynamics  equations  has  been 
obtained,  the  same  level  of  approximation  must  be  maintained  consistently 
throughout  any  subsequent  analysis,  in  the  course  of  such  analysis.  It  Is 
often  necessary  for  the  dependent  variable  In  a  term  to  be  replaced  with  an 
equivalent  expression.  It  Is  therefore  useful  to  note  that  the  dependent 
variables  In  second-order  terms  may  be  replaced  using  first-order 
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relations  without  decreasing  the  overall  level  of  the  approximation.  For 
example,  suppose  we  decided  to  recast  the  p'du/dt  term  In  Eq.  (B.1 ).  By 
using  Eq.  (B.6)  we  could  substitute  VP'/Pq  for  du/dt;  thus  an  equivalent 
second-order  expression,  (P7pq)  VP',  Is  formed.  If  we  had  substituted 
using  an  expression  from  Eq.  (B.1 1 ),  however,  a  higher-order  term, 
P'^VPq/Po,  would  be  Introduced.  This  term  would  have  to  be  dropped  In 
order  to  maintain  the  same  level  of  approximation.  This  extra  step  Is 
avoided  by  using  the  following  rule:  When  substituting  for  the  dependent 
variable  In  a  second-order  term,  use  a  first-order  relation  (see,  for 
example,  Lighthi  1 1 1 956). 

0.  Epilogue 

In  this  chapter  a  ranking  system  was  defined,  and  then  used  to 
simplify  the  lossless  hydrodynamics  equations  into  specialized  forms. 
These  forms  can  be  readily  combined  Into  ( I )  the  small-signal  wave 
equation  for  homogeneous  fluids,  (2)  the  small-signal  wave  equation  for 
inhomogeneous  fluids,  and  (3)  the  nonlinear  wave  equation  for 
inhomogeneous  fluids. 


CHAPTER  3 


LINEAR  GEOMETRICAL  ACOUSTICS 


A  Introduction 

In  this  chapter  linear  gec^etrical  acoustics,  commonly  referred 
to  as  ray  theory,  is  discussed,  and  the  mathematical  techniques  required 
for  nonlinear  geometrical  acoustics  are  introduced.  The  linear 
hydrodynamics  equations  for  a  lossless  inhomogeneous  medium  are 
combined  to  form  the  corresponding  wave  equation.  A  Galilean 
transformation  in  which  the  speed  of  the  moving  coordinate  system  is  the 
small-signal  sound  speed  Cq  is  introduced.  At  the  same  time  we  introduce 
what  is  called  the  geometrical  acoustics  assumption.  Use  of  the  Galilean 
transformation  and  geometrical  acoustics  assumption  enables  us  to 
separate  the  wave  equation  into  two  equations:  the  transport  equation  and 
the  elkonal  equation.  The  development  Is  performed  solely  In  the  time 
domain.  The  gradient  of  the  elkonal  is  related  to  the  variation  of  the  sound 
speed  in  the  medium.  For  the  case  of  a  stratified  ocean,  we  reduce  the 
elkonal  equation  to  an  expression  for  the  radius  of  curvature  of  the  ray 
path  The  transport  equation  Is  then  reduced,  via  two  transformations,  to 
the  form  of  the  first-order  plane  wave  equation  for  a  homogeneous  fluid. 
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The  analysis  starts  with  the  linear  hydrodynamics  equations 
which  include  inhomogeneity  terms.  By  the  end  of  the  derivation,  however, 
all  the  inhomogeneity  terms  have  been  neglected.  The  inhomogeneity  of  the 
medium  thus  enters  the  wave  equation  solely  through  the  dependence  of  the 
sound  speed  on  position.  Hence  the  results  of  this  chapter  are,  within  the 
ranking  scheme  of  Chapter  2,  valid  to  a  first-order  approximation. 


B.  Linear  Lossless  Wave  Eouation  for  Inhomogeneous  Fluids 

In  this  section  the  linear  lossless  hydrodynamics  equations  for  an 
inhomogeneous  fluid  are  combined  to  form  a  wave  equation.  First  the 
continuity  and  state  equations,  Eqs.  (2-B.  1 0)  and  (2-B.  1 3),  are  combined  to 
eliminate  the  dpVdt'  and  u-Vp^  terms.  The  time  derivative  of  the  result  is 

d^P'  du  0  ^  ^ 


(B.I) 


Use  of  the  momentum  equation,  Eq.  (2-B.  1 1 ),  eliminates  the  factor  du/dt 
term  in  Eq.  (B.  I ).  The  result  may  be  arranged  in  a  convenient  form 


v2p.  -  -L  - 

cl  at" 


1  d^p.  VPqVP*  VPq(VP*-c^Vp*) 


Po 


PqC? 


P'lVPol 


O'  +  £-  72p  _  ^  .  yp 

Po  0  ”  ° 


2P‘ 

Po 


(B.2) 


Since  the  three  terms  on  the  right-hand  side  of  Eq.  (B.2)  are  higher-order 
terms,  they  must  be  discarded  in  order  to  maintain  the  same  level  of 


approximation.  The  fourth  term  on  the  left-hand  side  of  Eg.  (B.2)  appears  to 
be  of  second-order  of  smallness,  an  inhomogeneity  term.  Recalling  the  rule 
mentioned  in  Chapter  2,  we  may  employ  a  first-order  relation  to  substitute 
for  the  dependent  variable  P'.  The  appropiate  first-order  relation  is  the 
linear  state  equation  for  a  homogeneous  fluid,  Eg.  (2-B.9).  Its  use  shows 
that  the  fourth  term  Is  actually  of  higher  order  and  must  therefore  be 
dropped.  The  effect  of  discarding  the  fourth  term  is  tantamount  to  having 
derived  the  wave  equation  disregarding  the  gravity  term  in  the  momentum 
equation,  Eg.  (2-A.3)  (Bergmann  1946).  By  a  more  physical  discussion, 
Bergmann  concluded  that  the  gravity  term  may  be  neglected  in  all  but  a  few 
extreme  situations,  such  as  the  propagation  of  very  low  frequency  signals. 

If  all  of  the  terms  mentioned  above  are  removed,  Eg.  (B.2)  may  be 
written  as  follows: 


-2-.  I  a^p'  ypo-vp' 

cj  dt^  Po 


*  0 


(B.3) 


The  first  and  second  terms  of  Eg.  (B.3)  constitute  the  classical  wave 
equation. 

The  third  term  is  present  because  of  the  inhomogeneity  of  the  fluid; 
this  term  was  also  discussed  by  Bergmann  (1946).  He  concluded  that  it  may 
be  neglected  if  the  density  gradient  is  sufficiently  small,  or  if  the 
frequency  of  the  signal  is  sufficiently  high.  (This  may  be  shown  by 
considering  an  example  situation  similar  to  that  used  below  in  Section  E.) 

In  fact  Brekhovskikh  and  Lysanov  (1982)  state  that  the  density  gradient 


term  may  be  neglected  for  frequencies  as  low  as  1  Hz.  If  the  density 
gradient  term  is  neglected,  Eq.  (B.3)  becomes 

V^P'  -  4  ^  “  0  ’ 

Co 

which  has  the  form  of  the  classical  wave  equation.  Thus  to  a  first-order 
approximation,  all  the  effects  of  the  fluid  inhomogeneity  enter  through  the 
variation  of  the  small-signal  sound  speed  Cq. 


The  assumption  that  we  are  about  to  make  is  referred  to  as  the 
geometrical  acoustics  assumption  (see,  for  example.  Landau  and  Llfshitz 
1959,  p.  256).  The  assumptions  made  up  until  now  can  be  seen  clearly  with 
the  aid  of  the  ranking  scheme;  only  first-order  terms  and  second-order 
inhomogeneity  terms  have  been  retained,  and  the  latter  have  been  found 
insignificant.  For  geometrical  acoustics  another  assumption  Is  needed.  We 
assume  that  the  surface  defined  by  an  arbitrary  wavefront  is  made  up  of 
many  small  segments  of  area,  each  of  which  may  be  regarded  as  plane.  This 
Is  a  reasonable  assumption,  since  we  are  concerned  with  long  range 
propagation.  At  long  ranges  the  curvature  of  a  spreading  wave  is  very 
small.  The  signal  associated  with  each  small  segment  of  area  follows  a 
path  that  is  perpendicular  to  that  segment.  The  path  is  called  a  ray;  see 
Fig.  3.1.  Within  the  geometrical  acoustic  assumption  the  travel  time 
associated  with  propagation  along  each  ray  is  defined  as  the  integral  of  the 
reciprocal  of  the  local  sound  speed  along  the  ray  path.  All  of  the  small 
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plane  surface  elements  that  make  up  the  wavefront  are  assumed  to  have  the 
same  phase  (travel  time),  and  hence  the  expression  equiphase  wavefronts. 

A  small  bundle  of  rays  may  be  thought  of  as  forming  a  tube.  Since  within 
the  geometrical  acoustics  assumption  no  propagation  takes  place  across 
rays,  a  propagating  disturbance  Is  bounded  by  the  walls  of  Its  ray  tube. 

The  geometrical  acoustics  assumption  has,  of  course,  a  range  of 
validity  (see,  for  example,  Frank,  Bergmann,  and  Yaspan  1969,  p.  63).  When  a 
ray  turns  appreciably  over  the  distance  of  a  wavelength,  the  equiphase 
wavefront  assumption  breaks  down  and  diffraction  has  a  strong  effect  on 
the  propagation  (signals  do  not  remain  confined  to  their  ray  tubes).  It  turns 
out  that,  for  a  ray  to  turn  appreciably,  the  sound  speed  must  vary  rapidly 
over  a  wavelength.  Hence,  for  geometrical  acoustics  to  be  a  valid 
assumption,  the  sound  speed  must  vary  slowly  on  a  wavelength  scale. 

Before  Imposing  the  geometrical  acoustics  assumption,  we  need  to 
Introduce  the  concept  of  retarded  time.  Consider  a  Galilean  transformation 
In  which  the  speed  of  the  moving  coordinate  system  Is  the  small-signal 
sound  speed  Cq.  In  this  system  the  new  time  t'  Is  termed  retarded  time. 

For  a  plane  wave,  retarded  time  Is  defined  as  follows: 

t'»t-x/CQ  ,  (CD 

where  x  Is  the  distance  the  wave  Is  propagated.  Note  that  x/Cq  Is  the  travel 
time  for  all  points  on  the  wavefront;  thus  plane  waves  are  equiphase 
wavefronts.  For  a  spherical  wave  retarded  time  Is 
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where  r  Is  the  radial  range  and  rp  Is  the  reference  range.  In  this  case  the 
travel  time  is  (r-roVc^.  Note  that  geometrical  acoustics  is  exact  in  the 
case  of  plane  and  spherical  waves. 

The  use  of  retarded  time  as  a  new  independent  variable  greatly 
simplifies  wave  equations.  To  illustrate  this  we  use  the  retarded  time 
variable  in  the  first-order  plane  wave  equation  for  a  wave  moving 
outwards  through  a  homogeneous  fluid.  A  first-order  plane  wave  equation 
may  be  formed  by  separating  the  commutable  derivatives  in  the  classical 
wave  equation  as  follows  (see,  for  example.  Pierce  1981 ,  p.  20); 

[d/di  -  Cq-'  d/axl  [d/di  *  Cq"’  d/dxj  P'  *  O  .  (C.3) 

The  first  term  in  square  brackets  may  be  thought  of  as  an  operator 
associated  with  Incoming  waves,  and  the  second  term,  with  outgoing  waves. 
Integration  with  respect  to  one  of  the  operators  leads  to  a  first-order 
wave  equation.  The  integration  constant  must  be  zero  to  satisfy  static 
conditions.  The  first-order  wave  equation  for  outgoing  waves  is  as 
follows; 

dP'/dt  ♦  Cq  dPVdx  -  0  .  (C.4) 


Substitution  of  the  outgoing  wave  function  P'(t  -  x/Cq)  into  Eq.  (C.4)  shows 
that  Eq.  (C.4)  is  the  correct  wave  equation.  If  new  Independent  variables  x, 
t',  where  t'  is  given  by  Eq.  (C.  1 ),  are  introduced.  Eq.  (C.4)  becomes 


Thts  Is  a  considerable  simplification.  The  physical  Interpretation  of 
Eq.  (C.5)  Is  simple:  Since  all  points  on  the  wave  move  with  the  same  velocity 
as  our  new  coordinate  system,  the  wave  appears  motionless.  Hence  the 
derivative  of  the  pressure  with  respect  to  x,  holding  t'  constant.  Is  zero. 

The  definition  of  retarded  time  In  an  Inhomogeneous  fluid  Is 
similar  to  Its  definition  for  plane  and  spherical  waves.  Because  the 
small-signal  sound  speed  may  change  along  the  ray  path,  the  travel  time  Is 
now  a  path  Integral: 

f  - 1  -  J  ds/c^  .  (C.6) 

where  ds  Is  an  Incremental  step  along  the  ray  path. 

We  now  state  a  general  definition  of  retarded  time  and  Introduce 
the  elkonal  t(r): 

t'«t-t(r)  .  (C.7) 

Note  that  the  elkonal  is  a  function  of  the  ray  path  position  vector  r;  It  Is 
not  a  function  of  t'.  The  elkonal  Is  defined  to  represent  the  travel  time,  or 
phase  shift,  associated  with  the  distance  between  the  origin  of  the  wave 
and  Its  current  position.  The  Integral  In  Eq.  (C.6)  plays  the  same  role.  It 
turns  out  that  In  the  case  of  Inhomogeneous  media,  this  definition  of  the 
elkonal  Is  equivalent  to  Imposing  the  geometrical  acoustics  assumption. 

Note  that  If  the  elkonal  t(r)  is  constant,  an  equiphase  wavefront  Is  defined. 
(For  this  reason,  equiphase  wavefronts  are  sometimes  referred  to  as 
elkonal  wavefronts.)  Equiphase  wavefronts  are  part  of  the  geometrical 


acoustics  assumption.  Thus  use  of  the  elkonal  In  our  definition  of  retarded 
time  means  that  the  geometrical  acoustics  assumption  is  invoked. 

The  rate  of  change  of  the  phase  of  the  equiphase  wavefront  may  be 
defined,  by  analogy  with  Eq.  (C.6),  as  follows: 


|Vt|  =  1/c, 


(C.8) 


This  definition  is  important  to  our  work  in  both  this  chapter  and  the  next. 

We  now  develop  the  relations  necessary  to  apply  the  Galilean 
transformation  based  on  Eq.  (C.7)  and  thus  to  introduce  the  geometrical 
acoustics  assumption.  First  the  effect  of  the  transformation  on  the  spatial 
and  temporal  derivatives,  V  and  d/at,  is  determined.  Careful  application  of 
the  chain  rule  leads  to 


dP'/at -» dP'/af  , 

VP'  -*  VP'  -  V tap* /at'  , 
v-u  -» v-u  -  vt-au/at'  , 


(C.9) 

(CIO) 

(CIO 


2  I  I 

V^P'-^|Vtp^  -  2^.(Vt'VP')  -  ♦V^P'  (C.12) 

at'^ 


As  is  shown  below,  the  terms  VtaP'/dt'  and  VP'  tnEq.  (CIO)  are, 
according  to  the  ranking  system  of  Chapter  2,  associated  with  first-order 
effects  and  second-order  effects,  respectively.  It  is  relatively  simple  to 
show  that  VtaP'/at'  is  associated  with  first-order  effects.  Recalling 
Eq.  (C.5)  we  note  that,  after  making  the  Galilean  transformation,  a  linear 
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progressive  plane  wave  appears  motionless  to  the  observer.  The  variation 
of  the  wave  occurs  solely  with  respect  to  the  retarded  time  variable  t'. 
Hence  the  derivative  of  a  field  variable,  such  as  dP'/dt',  is  a  first-order 
term. 

To  show  why  the  spatial  variation  VP'  is  associated  with 
second-order  effects,  we  again  recall  Eq.  (C.5).  In  that  equation  we  see  that 
the  spatial  variation  is  zero  in  the  case  of  a  small-signal  progressive  plane 
wave,  in  the  moving  coordinate  system  the  spatial  variation  term  can  have 
a  non-zero  value  only  if  a  wave  changes  its  shape  as  it  propagates.  Three 
causes  of  change  of  shape  are  the  decrease  in  amplitude  due  to  geometric 
spreading,  the  effects  of  inhomogeneity,  and  finite  amplitude  effects.  We 
are  interested  in  propagation  to  long  ranges,  well  into  the  farfield;  hence 
the  change  due  to  geometrical  spreading  is  small.  In  the  ranking  scheme  of 
Chapter  2,  the  terms  associated  with  small  effects  such  as  nonlinear 
distortion  are  ranked  as  second-order.  Accordingly  terms  containing  VP' 
are  ranked  as  second-order.  Note  that  we  assume  that  the  vector  VP'  is 
tangent  to  the  ray  path.  This  is  in  keeping  with  the  geometrical  acoustics 
assumption  wherein  the  propagation  of  the  signal  down  each  ray  tube  is 
treated  individually. 

Because  in  the  r,  t'  system  VP'  is  a  second-order  term,  the 
first-order  equivalent  of  Eq.  (C.10)  is 

VP’ = -VtdP'/dt’  .  (C.13) 

This  relation  is  useful  in  reducing  second-order  terms.  Recall  the 
substitution  rule  that  was  mentioned  in  Chapter  2.  For  example,  the 
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Galilean  transformation  of  the  second-order  quantity  Vp^-VP'  in  Eg.  (B.3) 
has  several  parts,  but  on  application  of  Eg.  (C.  1 3)  reduces  to 
(Vt-Vpo)dP7dt'. 

Equation  (C.13)  is  also  useful  for  maintaining  the  same  order  of 
approximation  in  first-order  relations.  For  example,  consider  the 
linearized  momentum  equation  for  a  homogeneous  medium.  Eg.  (2-B.6).  The 
Galilean  transformation  of  it  with  application  of  Eg.  (C.  1 3)  in  place  of 
Eg.  (CIO)  Isas  follows: 

Podu/df  =  VtdP'/dt'  .  (C.14) 

Because  Is  not  a  function  of  t‘,  we  may  easily  Integrate  Eg.  (C.  14)  with 
respect  to  tV  Note  that  the  Integration  constant  must  be  zero  in  order  to 
satisfy  static  conditions.  The  result  is 

u»(P7po)Vt  .  (C.I5) 

This  expression  may  be  thought  of  as  a  generalization  of  the  progressive 
plane  wave  impedance  relation,  P'  «  PoCqU.  Equations  (C.  1 3)  and  (C.  1 5)  are 
used  In  the  next  chapter  to  simplify  nonlinear  terms. 

We  now  apply  the  Galilean  transformation  to  the  small-signal 
wave  equation  for  inhomogeneous  media.  Eg.  (B.4).  Using  the  transforms  of 
the  spatial  and  temporal  derivatives,  Egs.  (C.9)  and  (C.  12).  we  obtain  the 
following: 


I 

( 

( 

t 
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=  V^P'.  (C.16) 

in  order  to  assess  the  relative  Importance  of  the  term  on  the 
right-hand  side  of  Eq.  (C.16),  V^P',  we  need  to  understand  its  physical 
significance.  The  Laplacian  of  the  acoustic  pressure  P'  is  defined  to  be  the 
divergence  of  the  gradient  of  P';  that  is,  7^P'  as  V  VP'.  It  has  already  been 
argued  that,  in  the  moving  coordinate  system,  VP'  is  related  to 
second-order  effects,  the  variation  of  the  shape  of  propagating  wave.  Hence 
V*  VP'  is  the  rate  of  change  of  the  variation,  in  the  case  at  hand, 
small-signal  propagation,  the  variation  of  the  propagating  wave  is  due 
solely  to  geometric  spreading.  In  geometrical  acoustics  the  ray  paths  are 
not  permitted  to  turn  appreciably  over  a  wavelength;  hence  the  spreading 
takes  place  gradually.  The  rate  of  change  of  the  spreading  is  therefore 
assumed  to  be  zero.  In  the  case  of  finite  amplitude  waves,  the  situation  is 
not  so  simple.  As  already  noted  the  term  VP'  is  assumed  to  be  associated 
with  the  second-order  effects  of  finite  amplitude  and  inhomogenelty;  it  Is 
also  assumed  to  be  tangent  to  the  ray  path.  We  now  assume  the  distortion 
effects  occur  slowly  and  without  sudden  changes.  With  this  assumption  the 
term  V^P'  can  be  neglected  since  |  VP'  |  »  IV^P'  |.  Note  that  in  other 
situations,  such  as  those  in  which  the  rays  turn  rapidly,  the  term  V^P'  is 
important  since  it  accounts  for  diffraction  (Pelinovsky,  Petukov,  and 
Fridman  1979). 

If  the  right-hand  member  is  neglected,  Eq.  (C.  16)  may  be  expressed 
as  follows: 


(C.t7) 


d_ 

df 


=  0 


We  now  split  Eq.  (C.  1 7)  into  two  parts.  To  do  this  we  use  Eq.  (C.8)  which  was 
obtained  from  the  geometrical  acoustics  assumption,  and  which  shows  that 
the  first  part  of  Eq.  (C.  1 7)  vanishes, 

l/q,2- |vir|2*o  .  (C.18) 


Equation  (C.18)  is  referred  to  as  the  eikonal  equation.  As  is  shown  in  the 
next  section,  the  eikonal  equation  dictates  the  acoustic  ray  paths.  The 
remainder  of  Eq.  (C.  1 7)  may  be  integrated  once  with  respect  to  the  retarded 
time  variable  t'.  Noting  that  the  integration  constant  must  be  zero  in  order 
to  satisfy  static  conditions,  we  arrive  at  the  following: 

V*(VtP’2)  =  0  .  (C.19) 


Equation  (C.19)  is  called  the  transport  equation  and  is  equivalent  to 
Eq.  (8-5.13)  in  the  book  by  Pierce  (1981).  It  turns  out  that  this  equation 
governs  the  variation  of  the  acoustic  pressure  along  the  ray  path. 


D.  Ray  Paths  from  the  Eikonal 

In  this  section  it  is  shown  that  the  eikonal  equation,  Eq.  (C.  18), 
defines  the  acoustic  ray  paths.  (A  ray  path  and  the  coordinate  system  are 
shown  in  Fig.  3.2.)  First  the  connnection  between  the  eikonal  and  the  ray 
path  is  found.  Next  the  ray  path  is  linked  to  the  spatial  variation  of  the 
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sound  speed.  We  then  combine  the  two  ideas  and  find  a  general  relation  for 
the  ray  paths  in  terms  of  the  gradient  of  the  eikonal,  the  gradient  of  the 
sound  speed,  and  the  tangent  to  the  ray  path.  The  general  relation  is  then 
simplified  for  the  case  of  a  stratified  medium. 

A  connection  between  the  eikonal  and  the  ray  path  is  sought.  The 
eikonal  defines  an  equiphase  surface;  thus  Vi'  is  by  definition,  normal  to 
that  surface  (see,  for  example,  Thomas  1968,  Sec.  15.6).  Since  Eg.  (C.18) 
gives  the  square  of  the  magnitude  of  Vi',  it  is  easily  seen  that 


CqV^  “  n  , 


(D.l) 


where  n  is  the  unit  normal  to  the  surface.  It  is  known  from  vector  calculus 


that  the  unit  tangent  to  the  ray  path  T  is  defined  as  the  derivative  of  the 
ray  path  position  vector  r  with  respect  to  the  ray  path  lengths. 


T  s  dr/ds  . 


(D.2) 


Since  the  ray  path  is  normal  to  the  wavefront,  the  unit  tangent  to  the  ray 
path  and  the  unit  normal  to  the  wavefront  are  equal, 


T  =  n 


(D.3) 


(Note  that  they  are  not  necessarily  equal  in  the  case  of  a  moving  medium; 
see,  for  example,  Pierce  1 98 1 ,  p.  37 1 .)  Combining  Eqs.  (D.  1 )  and  (D.2)  yields 
the  following  relation: 


CqV^  -  dr/ds  . 


(D.4) 


•m-Cv-C 
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The  gradient  of  the  eikonal  is  thus  closely  related  to  the  tangent  to  the  ray 
path. 

An  equation  that  links  the  ray  path  with  the  variation  of  the  sound 
speed  in  the  medium  is  now  sought.  From  the  Frenet-Serret  formulas  (see, 
for  example.  Sokolnikoff  and  Redheffer  1958,  Sec.  4.1  i ),  we  know  that  a  ray 
path  may  be  thought  of  as  a  sequence  of  arcs  of  radius  R,  wherein  R  varies 
as  a  function  of  position  along  the  ray  path.  The  following  equations  are 
part  of  the  Frenet-Serret  formulas: 


KN  =  dT/ds  . 

(D.5) 

Rsl/K 

(D.6) 

where  K  is  the  curvature  and  N  is  the  principal  normal.  The  principal 
normal  Is  defined  to  be  normal  to  T  and  point  towards  the  center  of 
curvature.  From  the  above  relations  the  radius  of  curvature  R  may  be 
expressed  as  follows: 


R*|dT/ds|*’  .  (D.7) 

Because  of  its  close  relation  to  the  ray  path,  let  us  find  an 
expression  for  R  or,  equivalently,  an  expression  for  dT/ds.  Use  of  Eq.  (D.2) 
yields 


dT/ds  *  d(dr/ds)/ds  .  (D.8) 

Noting  Eq.  (D.4)  and  the  relation  d/ds  -  (dr/ds  •  V),  we  see  that  Eq.  (D.8) 
becomes 
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dl/ds  =  (dr/ds  •  VKCoVt)  .  (D.9) 

Reappllcatlon  of  Eq.  (D.4)  gives 

dJ/ds  =  (c^Vt  •  V)  (CoVi^)  .  (D.  1 0) 

Use  of  the  chain  rule  and  the  following  relation  (see,  for  example, 
Gradshteyn  and  Ryzhik  1 980,  Eq.  1 0.3 1 .3), 


inEq.  (D.IO)  leads  to 

dT/dS  =  Vt(CoVt- VCq)  +  (Co2/2)  V(Vt*  Vt)  .  (DID 

Use  of  the  eikonal  equation,  Eq.  (C.T8),  and  the  chain  rule  in  Eq.  (D.  1 1 }  yields 

dT/ds  =  CoVt(  Vt-Vc^)  -  (  Vc^>)/Cq  .  (D.  1 2) 

Recalling  Eq.  (D.7)  we  note  that  Eq.  (D.I2)  connects  the  radius  of  curvature  R 
to  the  variation  in  the  sound  speed  VCq.  Since  the  ray  paths  are  arcs  of 
radius  R,  Eq.  (0. 12)  may  be  thought  of  as  connecting  the  ray  paths  with  the 
variation  of  the  sound  speed.  Equation  (D.  1 2)  is  equivalent  to  Eq.  (66.6)  in 
the  book  by  Landau  and  Lifshitz  ( 1 959). 

A  physical  understanding  of  Eq.  (D.  12)  may  be  obtained  by 
considering  a  simple  example.  Let  the  medium  be  azimuthal  1y  symmetric 
(thus  restricting  ourselves  to  the  x,z  plane)  and  stratified,  that  is,  the 
sound  speed  profile  does  not  vary  in  range.  (Recall  that  k  is  in  the  direction 


of  increasing  depth.)  An  expression  for  the  radius  of  curvature  R  can  be 
found  by  reexamining  Eq.  (D.  1 2).  Under  conditions  stated  above,  the  terms  in 
Eg.  (D.  12)  may  be  defined  as  follows: 

VCo  =  kh  ,  (D.13) 

dr/ds  -  i  cos  0  +  k  sin  e  ,  (D.  1 4) 

where  h  is  the  gradient  of  the  sound  speed  profile  dCQ/dz,  6  Is  the  acute 
angle  between  the  ray  path  and  the  horizontal  (see  Fig.  3.2),  and  I  and  k  are 
unit  vectors  in  the  Cartesian  system.  If  Eq.  (D.4)  is  used,  Eq.  (D.  1 4)  may  be 
recast  in  the  following  form: 

Vf  =  I  Cq'>  cos  0  +  k  Cq*’  sin  0  .  (D.15) 

Substituting  Eq.  (D.  1 3)  and  (D.  1 5)  into  (D.  1 2)  yields  the  following: 

aT/ds  =  h  Cq~’  cos  0  (I  sin  0  -  k  cos  0)  .  (D.16) 

Let  us  interpret  Eq.  (0.16).  From  a  geometrical  consideration  of 
the  equation  we  see  that,  if  the  gradient  h  Is  positive,  the  center  of 
curvature  is  above  the  current  ray  path  position.  The  converse  is  true 
when  h  is  negative.  Thus  the  acoustic  rays  bend  toward  the  region  of  the 
slower  sound  speed.  From  the  reciprocal  of  the  magnitude  of  Eq.  (D.  1 6),  we 
obtain  the  radius  of  curvature: 


The  quantity  Cq/cos  0  In  Eq.  (0.17)  is,  according  to  SnelVs  law,  a  constant. 
From  Eq.  (0.17)  and  Snell's  law  we  see  that,  if  h  is  constant,  the  radius  of 
curvature  R  must  be  constant.  Therefore  rays  in  a  constant  gradient 
medium  are  circular  arcs.  Equation  (0.17)  appears  in  a  variety  of 
references  (see,  for  example.  Officer  1958,  Eq.  2-82). 

In  summary  it  has  been  shown  that,  from  the  eikonal  equation  and 
the  relations  of  vector  calculus,  a  general  expression  for  the  ray  paths  may 
be  found.  Moreover,  in  the  case  of  a  constant  gradient  medium,  the  ray 
paths  are  circular  arcs. 

E.  Acoustic  Pressure  from  the  Transport  Equation 

In  this  section  it  is  shown  that  the  transport  equation,  Eq.  (C.19), 
may  be  placed  in  the  form  of  the  first-order  plane  wave  equation  in  the 
moving  coordinate  system,  Eq.  (C.5).  It  turns  out  that  the  transport 
equation  governs  the  acoustic  pressure  at  any  point  along  the  ray  path,  and 
that  the  acoustic  pressure  is  related  to  the  area  of  the  ray  tube.  To  show 
this  we  introduce  the  concept  of  ray  coordinates  and  use  an  expression  that 
is  required  in  our  development  of  nonlinear  geometrical  acoustics.  (The 
details  of  the  transformation  from  Cartesian  coordinates  to  ray 
coordinates  are  included  in  Appendix  A.)  The  three  ray  coordinates  are  the 
ray  path  length  s  and  the  initial  ray  launch  angles  4>  and  i|i;  see  Fig.  (3.2). 

To  reduce  the  transport  equation  we  need  an  expression  for 
The  following  relationship  is  derived  in  Appendix  A: 
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JLA/^\ 
^0  1^0/ 


(E.1) 


where  Aq  Is  the  ray  tube  area.  It  Is  shown  below  that  the  dependence  of 
on  the  variation  in  the  sound  speed  may  be  neglected.  We  first 
rearrange  Eg.  (E.1 )  as  follows. 


y2^ 


.  _Lf^!Og  A(^\  -  £oi  A/5o_\] 

Co[Aq  dS\AQ,/  Cq  dsyCog/J 


(E.2) 


where  Aq,  and  Cq,  are  the  ray  ti^  area  and  the  sound  speed  at  the 
reference  ray  path  length  ( I  meter  in  this  example).  Consider  vertical 
propagation  in  a  stratified  ocean  5(X)0  m  deep.  Let  the  sound  speed 
increase  linearly  with  depth  from  1400  m/s  to  Cq,  *  1600  m/s;  i.e., 

Cq  *  1400  ♦  0.04  z,  where  z  is  the  depth.  This  is  a  larger  variation  than  is 
usually  expected  in  the  ocean  and  may  be  regarded  roughly  as  an  upper 
bound.  In  the  case  of  vertical  propagation,  d/ds  ■  d/dz  and  the  ray  tube 
areas  may  be  approximated  by  using  spherical  spreading.  The  ratio  of  the 
ray  tube  areas  Aq/Aq,  is  z^,  hence  d(Ao/Ao,)/dz  =  2z.  If  the  propagation 
path  is  the  full  ocean  depth,  5000  m,  the  term  is  equal 

to  4x  10“®^.  Using  the  conditions  given  above,  the  term  (Co5/Co)d(Co/Co8)/ds 
can  be  seen  to  be  approximately  3  x  lO'^.  Therefore,  even  when  the 
variation  of  the  sound  speed  is  unusually  large,  the  term 
(Co,/Co)d(Co/Co5)/ds  is  small  with  respect  to  (Aos/AQ)d(Ao/Aog)/ds.  The 
variation  in  the  sound  speed  is  therefore  neglected  in  Eg.  (E.2),  and  Eg.  (E.1 ) 
reduces  to 


-  (CqAo)'’  dV^s 


(E.3) 


One  other  relation  Is  required  to  reduce  the  transport  equation. 
Starting  with  the  expression  d/ds  « (dr/ds  *  V),  and  recalling  Eq.  (D.4),  we 
see  that 


Vi'*  V()  =  (l/Co)d()/ds  .  (E.4) 

Stmpliricatlon  of  the  transport  equation.  Eq.  (C.I9).  is  obtained  by 
forcing  It  toward  the  form  of  Eq.  (C.S).  Substitution  of  Eqs.  (E.3)  and  (E.4)  In 
the  transport  equation  gives  the  following: 


J!L  .  0 

ds  2Ao  ds 


(E.5) 


Equation  (E.5)  may  be  expressed  as  a  perfect  differential  by  multiplying 
both  sides  by  and  rearranging, 

d(P'Ao’'2)/as  =  0  .  (E.6) 

A  physical  Interpretation  of  Eq.  (E.6)  is  that  P'Aq'^^  is  constant  along  the 
ray  tube,  a  result  consistent  with  the  assumption  that  the  energy  Is 
constant  along  the  tube.  The  fact  that  P'Aq'^^  Is  constant  suggests  a  new 
dependent  variable,  one  for  which  the  amplitude  Is  not  affected  by 
geometric  spreading.  To  make  the  new  dependent  variable  have  the  same 
units  as  P',  It  Is  defined  as  follows: 


w  S  (  VAqs)’'^  P' 
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If  this  transformation  is  used,  Eq.  (E.6)  becomes 

dW/ds  -  0  ,  (E.8) 

which  has  the  desired  form  of  Eq.  (C.5) 

We  now  obtain  the  solution  of  the  problem  of  propagation  of 
small-signal  waves  through  an  inhomogeneous  medium.  The  analogous 
solution  for  a  plane  wave  in  a  homogeneous  fluid  is  as  follows.  Given  the 
boundary  condition 

P' =g(t)  =  g(t’)  at  x  =  0  .  (E.9) 

the  solution  is 

P’=g(t’)  ,  (E.IO) 

where  t'  is  given  by  Eq.  (C.  1 ).  For  the  problem  of  waves  propagating  in  an 
inhomogeneous  medium,  the  boundary  condition  is 

P' =g(t)  =  g(t')  at  s  =  So  .  (Eli) 

Although  the  boundary  condition  does  not  have  quite  the  same  form  as 
Eq.  (E.9),  the  form  can  be  made  exactly  the  same  by  defining  a  new 
independent  variable.  The  new  variable  should  not,  however,  alter  the  form 
of  Eq.  (E.8).  The  desired  transform  is 

Z  =  s-So  .  (E.12) 


(Note  that  Z  is  not  the  same  as  z,  the  depth.)  The  boundary  condition  is  now 
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P'=g(t)  =  9(t')  at  Z  =  0  ,  (E.13) 

where 

t'  =  t-z/Co 

Equation  (E.8)  now  becomes 

dW/dZ  =  0  .  (E.15) 

Equation  (E.  1 5)  is  in  the  form  of  Eq.  (C.5).  Since  Eq.  (E.  1 3)  is  equivalent  to 
Eq.  (E.9),  the  solution  of  the  problem  of  propagation  in  an  inhomogeneous 
medium  is,  by  analogy, 

W  =  g(t‘)  .  (E.16) 


where  Z  and  W  are  defined  in  Eqs.  (E.  12)  and  (E.7). 

Thus  the  linear  transport  equation  for  an  Inhomogeneous  fluid, 
Eq.  (C.19),  has  been  transformed  into  an  equation  in  the  form  of  the 
first-order  plane  wave  equation  for  a  homogeneous  medium  by  making  two 
simple  transformations:  one  on  the  dependent  variable  P',  Eq.  (E.7),  and  the 
other  on  the  independent  variable  s,  Eq.  (E.  12).  The  main  reason  for 


introducing  the  W  and  Z  transformations  is  to  prepare  the  reader  for 
similar,  but  more  complicated,  transforms  in  Chapter  4. 


CHAPTER  4 


NONLINEAR  GEOMETRICAL  ACOUSTICS 


This  chapter  progresses  in  a  fashion  only  slightly  different  from 
that  of  the  previous  chapter.  In  Chapter  3  the  linear  lossless 
hydrodynamics  equations  for  an  inhomogeneous  fluid  were  combined  to  form 
the  corresponding  wave  equation.  Next  a  Galilean  transformation,  based  on 
the  use  of  the  retarded  time  t',  and  the  geometrical  acoustics  assumption 
were  introduced.  The  Galilean  transformation  was  then  applied  to  the  wave 
equation.  In  this  chapter  we  start  with  the  nonlinear  lossless 
hydrodynamics  equations.  They  are  not,  however,  combined  to  form  a 
nonlinear  wave  equation,  but  instead  only  partially  combined.  The  Galilean 
transformation  and  the  geometrical  acoustics  assumption  are  then 
introduced.  The  hydrodynamics  equations  are  transformed  and  then 
combined  to  form  a  wave  equation  in  the  moving  coordinate  system.  From 
this  point  on,  the  development  parallels  that  in  Chapter  3  After  the  V^p' 
term  is  dropped,  the  transformed  wave  equation  is  separated  into  the 
eikonal  equation  and  the  transport  equation.  The  equations  are  then 
discussed.  After  simplification  the  transport  equation  is  reduced  to  a 
first-order  (nonlinear)  wave  equation.  This  equation  has  the  same  form  as 
the  first-order  equation  for  plane  waves  of  finite  amplitude  in  a 


homogeneous  medium.  The  reason  we  depart  from  the  procedure  followed  In 
Chapter  3  is  that  It  is  simpler  to  do  so.  It  is  easier  to  transform  and 
approximate  the  hydrodynamics  equations  and  then  combine  them  to  obtain  a 
wave  equation  than  vice  versa. 

In  this  chapter  the  standard  geometrical  acoustics  assumption 
must  be  broadened,  in  particular  the  acoustic  ray  paths  are  assumed  to  be 
unaffected  by  the  signal;  that  is,  self-refraction  does  not  occur.  This 
assumption  means  that  acoustic  ray  path  equations  developed  in  the 
previous  chapter  may  be  used  even  for  signals  of  finite  amplitude. 
Specifically,  the  eikonal  equation,  Eq.  (3-C.  18),  and  the  equation  developed  in 
Appendix  A,  Eq.  (3-E.3),  are  still  valid. 

The  level  of  approximation  used  in  this  chapter  differs  from  that 
used  in  the  previous  chapter.  There  the  only  second-order  terms  that  were 
considered  were  the  inhomogeneity  terms.  In  the  end  even  these  terms 
were  found  to  be  so  small  that  they  were  neglected.  Within  the  ranking 
system  of  Chapter  2,  this  means  that  the  results  of  the  previous  chapter 
are  valid  to  a  first-order  approximation  only.  In  this  chapter  we  are 
specifically  interested  in  examining  the  effects  of  another  type  of 
second-order  term,  nonlinear  terms.  In  order  to  be  consistent  in  our  level 
of  approximation,  terms  due  to  both  inhomogeneity  and  nonlinearity  are 
included  throughout  the  entire  derivation.  Thus  within  the  ranking  system 
of  Chapter  2,  the  results  obtained  In  this  chapter  are  valid  to  a 
second-order  approximation. 

As  was  mentioned  in  Chapter  1,  no  new  material  Is  presented  in 
this  chapter.  The  following  derivation  closely  follows  that  of  Pelinovsky, 
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Petukhov,  and  Fridman  (1979),  whose  work  is  in  turn  related  to  that  of 
Ostrovsky,  Pelinovsky,  and  Fridman  (1976).  Both  of  the  cited  papers  are, 
however,  terse.  So  much  so  in  fact  that  it  was  decided  that  the  full 
derivation  of  nonlinear  geometrical  acoustics  should  be  presented  here. 
That  decision  in  turn  led  to  the  inclusion  of  much  background  material  that 


appears  in  Chapters  2  and  3. 


The  lossless  nonlinear  hydrodynamics  equations  for  an 
inhomogeneous  fluid  were  developed  in  Chapter  2.  To  combine  the  nonlinear 
continuity  and  momentum  equations,  Eqs.  (2-A.8)  and  (2-B.  1 ),  we  take  the 
time  derivative  of  continuity  equation,  the  divergence  of  the  momentum 
equation,  and  substitute  the  former  into  the  latter  via  the  V*Podu/dt  term. 
Placing  all  nonlinear  terms  on  the  right-hand  side,  we  obtain  the  following 


equation: 


dt^  Po 


-V  P’*— =  V- Po(u'V)u-^vP--^  (A.1) 


p’  d(p'u) 

—  VP’ - - — - 

Po  dt 


To  arrive  at  Eq.  (A.1),  we  neglected  several  third-order  terms;  terms 
involving  the  Laplacian  of  environmental  parameters  such  as  V^Pq  ,  and 
terms  involving  the  products  of  gradients  of  environmental  parameters 
such  as  VPo'VPq.  It  was  also  necessary  to  use  the  linear  momentum 
equation  for  a  homogeneous  fluid,  Eq.  (2-B.6),  to  simplify  one  of  the 


second-order  terms. 


The  next  step  Is  to  eliminate  p'  and  u  from  Eq.  (A.  I ).  To  develop 
some  of  the  relations  required  to  do  this,  we  turn  our  attention  to  the 
equation  of  state,  Eq.  (2- A.  1 ).  The  second-order  Taylor  series  expansion  of 
the  equation  of  state  Is  given  In  Eq.  (2-B.  15).  We  now  find  the  temporal 
derivative  of  Eq.  (2-B.  15)  and  the  gradient  of  Eq.  (2-A.  1 ): 


at  ‘'"at  «  lap  1 


t.X4>  -  pQ 


-  (u  V)Jf^]  -  (u  V)xM 

Wp.X.t-to  ^  Jp^.X-Xt 

VP-c2vp.vJ|ft  .  vx® 

^‘^^Jp.X4-to  ^  ^-'pA.X-Xo 


(A.2) 


(A3) 


In  arriving  at  Eq.  (A.2),  we  have  made  use  of  the  fact  that  the  entropy  and 
salinity  of  a  material  particle  remain  constant;  see  Eqs.  (2-B.2)  and  (2-B.4). 
We  have  also  used  the  definition  of  the  small-signal  sound  speed, 

Eq.  (2-B.8). 

To  reduce  the  effort  required  to  combine  Eqs.  (A.1 ),  (A.2),  and 
(A.3),  we  make  use  of  the  Galilean  transformation  and  geometrical 
acoustics  assumption,  which  are  discussed  at  length  in  the  previous 
chapter. 


B.  Geometrical  Acoustics  Assumption 


The  Galilean  transformation  and  the  geometrical  acoustics 


assumption  are  used  in  the  previous  chapter  to  reduce  the  linear  wave 
equation  to  the  elkonal  equation  and  a  first-order  wave  equation  that  has 
the  form  of  Eq.  (3-C.5).  In  this  chapter  we  transform  the  Ingredients  of  the 
wave  equation,  that  Is,  the  nonlinear  hydrodynamics  equations  for 
Inhomogenenous  fluids.  The  equations  are  then  combined  to  form  a  wave 
equation  In  the  moving  coordinate  system. 

Before  the  Galilean  transformation  Is  applied  to  the  nonlinear 
hydrodynamics  equations.  It  Is  useful  to  review  the  transformation  as 
applied  to  the  equation  for  finite  amplitude  plane  waves  propagating 
through  a  homogeneous  medium.  The  equation  Is  as  follows  (see,  for 
example,  Blackstock  1972): 

P'i*CoP'„-(p/PoC(,2)P'P\»0  ,  (B.1) 

where  <  he  subscripts  t  and  x  Indicate  derivatives  with  respect  to  time  and 
distance,  respectively,  and  p  Is  the  coefficient  of  nonlinearity 
(see  Appendix  B).  Equation  (B.  1 )  is  the  nonlinear  extension  of  Eq.  (3-C.4). 
Introduction  of  the  retarded  time  variable  for  plane  waves,  Eq.  (3-C.  1 ), 
transforms  Eq.  (B.  1 )  into  the  following: 

p',  -  <P/PoCo’>  P'P'i'  =  0  (B2) 

Equation  (B.2)  is  the  equation  for  outgoing  plane  waves  in  the  second 
approximation.  It  Is  the  nonlinear  extension  of  Eq.  (3-C.5).  Equation  (B.2) 
has  been  analyzed  extensively  in  the  past  and  its  solution,  called  the 
Eamshaw  solution,  is  well  known  (see,  for  example,  Blackstock  1972)  Use 
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of  the  Galilean  transformation  and  the  geometrical  acoustics  assumption 
enables  the  nonlinear  hydrodynamics  equations  to  be  put  Into  the  form  of 
Eg.  (B.2);  then  the  problem  of  propagation  In  an  Inhomogeneous  medium  will 
have  been  solved. 

The  Galilean  transformation  is  now  applied  to  Eqs.  (A.1),  (A.2),  and 
(A.3).  The  new  time  variable  t'  is  defined  by  Eq.  (3-C.7).  Expressions  for 
various  spatial  and  temporal  derivatives  are  given  In  Eqs.  (3-C.9)  through 
(3-C.12).  Since  the  Galilean  transformation  is  a  mathematical  coordinate 
transformation,  It  Is  unaffected  by  the  Inclusion  of  nonlinear  terms. 
Equation  (A.  1 )  becomes 


af’  '  '  at'2  if  ''  ^af  ''  p  at’ 


at’ 


at' 


at’ 


-  vt  •  Po(u- vt)—  ♦ 

...  at’ 


at' 


(B.3) 


To  arrive  at  Eq.  (B.3),  we  have  used  the  first-order  expression  for  the 
gradient,  Eq.  (3-C.13),  in  nonlinear  and  Inhomogeneity  terms.  The  spatial 
and  temporal  derivatives  of  the  equation  of  state,  Eqs.  (A.2)  and  (A.3) 
become,  respectively. 


(B.4) 


7P„  *  7P-  -  7*1^  -  c’  (vpo  ♦  Vp'  -  Vi-  ^  I 
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I  ^0 
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at'  /[axj 


(B.5) 


P.I.X  ■  Xq 


Equations  (B.4)  and  (B.5)  are  now  rearranged  In  anticipation  of 
combining  them  with  Eq.  (B.3).  First  consider  Eq.  (B.4);  take  Its  time  (f) 
derivative,  remembering  that  Vi'  Is  Independent  of  t'.  Use  of  the  linear 
momentum  equation  In  the  moving  coordinate  system,  Eq.  (3-C.15),  the 
elkonal  equation,  Eq.  (3-C.18),  and  the  linear  equation  of  state,  Eq.  (2-B.9), 
results  In  the  following  equation; 


I  d^P' 

P„cjafraf/ldxf 


'p.w*»o 


at'  \  at'  /  lau  . ,  . 

<Q  \  !  K  '>)  p.X.t»to 


(B.6) 


Multiplication  of  Eq.  (B.5)  by  Vi'(p'/Po)  places  It  In  a  more  readily  usable 
form.  The  same  relations  used  In  arriving  at  Eq.  (B.6).  namely.  Eqs.  (3-C.  14). 
(3-C.18).  and  (2-B.9)  are  then  used  to  simplify  the  result.  Rearranging  and 
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taking  the  derivative  with  respect  to  t'  yields 


ap'  vt-yPp  vpo'V't'ap'  i  a  /  dt'\M 
Po  ■  Po  «''poC^Pn 


p.t,X  “  Xq 


(B.7) 


Equations  (B.6)  and  (B.7)  may  be  combined  with  Eq.  (B.3).  First 
simplify  the  right-hand  side  of  Eq.  (B.3)  using  the  same  techniques  used  to 
arrive  at  Eqs.  (B.6)  and  (B.7).  Substitution  of  Eqs.  (B.6)  and  (B.7)  into  the 
left-hand  side  of  Eq.  (B.3)  leads  to 


_d_ 

dV 


VPo-V^\ 
Po  / 


P’ 


-  V^P' 
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j 

(B.8) 


Where  the  coefficent  of  nonlinearity  p  Is  defined  as  follows: 


(B.9) 


Equation  (B.8)  Is  similar  to  the  linear  geometrical  acoustics  equation. 

Eq.  (3-C.16).  The  main  difference  is  that  Eq.  (B.8)  has  a  nonlinear  term  on 
the  right-hand  side.  Another  difference  is  that  the  left-hand  side  contains 
an  Inhomogenelty  term,  namely,  (Vpo-Vt/po)P'.  The  ancestor  of  this  term 
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is  the  third  term  in  Eq.  (3-B.3),  there  referred  to  as  the  density  gradient 
term.  Neglected  in  Chapter  3  as  being  unimportant  except  for  exceedingly 
low  frequencies,  the  density  gradient  term  is  retained  here.  In  this  chapter 
we  are  specifically  interested  in  examining  finite  amplitude  effects.  In 
order  to  be  consistent  in  our  level  of  approximation,  we  maintain  the 
second-order  terms  due  to  both  inhomogeneity  and  finite  amplitude  effects 
throughout  the  entire  derivation. 

Finally,  as  in  linear  geometrical  acoustics,  we  discard  V^P'.  It 
was  noted  in  the  previous  chapter  that  this  term  is  associated  with  the 
rate  of  change  of  the  distortion.  Since  the  distortion  is  assumed  to  occur 
slowly,  the  rate  of  change  of  distortion  must  be  very  small  and  is  therefore 
neglected.  The  final  form  of  the  nonlinear  geometrical  acoustics  equation 
is  as  follows; 


C.  Reduction  of  the  Nonlinear  Geometrical  Acoustics  Equation 

The  techniques  used  in  this  section  to  reduce  the  nonlinear 
geometrical  acoustics  equation  are  the  same  as  those  used  in  Chapter  3  to 
reduce  the  linear  geometrical  acoustics  equation.  Since  we  assume  that  the 
geometrical  acoustics  assumption  is  valid  in  the  case  of  finite  amplitude 
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waves,  we  may  use  the  elkonal  equation,  Eq.  (3-C.18),  to  eliminate  the  first 
two  terms  on  the  left-hand  side  of  Eq.  (B.  10).  The  remaining  terms 
constitute  the  transport  equation  which  in  this  case  is  nonlinear.  The 
procedure  used  in  Chapter  3  to  reduce  the  transport  equation  to  a 
first-order  wave  equation  is  now  repeated.  In  this  case  the  first-order 
wave  equation  has  the  form  of  Eq.  (B.2).  Noting  that  is  independent  of 
t',  we  integrate  the  remaining  terms  in  Eq.  (B.IO)  once  with  respect  to  t'. 
The  integration  constant  must  be  zero  in  order  to  satisfy  static  conditions. 
The  result  is 


2Vt 


\  Po  /  at 


tc,i) 


Use  of  the  same  equations  used  to  reduce  the  linear  transport  equation, 
namely  Eqs.  (3-E.3)  and  (3-E.4)  and  the  eikonal  equation,  Eq.  (3-C.18),  enables 
Eq.  (C.  1 )  to  be  rewritten  as  follows; 


iE! 

ds 


2Aods\Co/  2pods 


(C.2) 


The  transport  equation,  Eq.  (C.2),  would  be  very  similar  to  the 
finite  amplitude  plane  wave  equation,  Eq.  (B.2),  if  the  middle  term  could  be 
forced  to  vanish.  To  force  Eq.  (C.2)  towards  the  form  of  Eq.  (B.2),  we 
introduce  a  transformation  of  the  dependent  variable; 
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Where  k  is  a  function  (not  a  constant)  to  be  defined  in  such  a  way  that  the 
boundary  condition  at  the  source  transforms  conveniently,  that  is,  W  =  P*  at 
the  source.  Use  of  Eq.  (C.3)  in  Eq.  (C.2)  yields 


Setting  the  middle  term  to  zero  yields  a  differential  equation  in  k.  which 
may  be  cast  in  the  form  of  a  perfect  differential 


a(£n  (AQ/pQCQk^)]/ds  =  0 


(C.5) 


The  solution  of  Eq.  (C.5)  is 


k  (QqCq/Aq)^^^  =  constant 


(C.6) 


We  now  choose  the  integration  constant  to  suit  the  boundary  conditions  at 
the  source.  Accordingly  the  definition  of  k  is  as  follows: 


k  -  (AqPqjCqj/AqjPqCq)’ 


(C.7) 


where  the  density,  the  small-signal  sound  speed,  and  the  ray  tube  area  at 
the  source  position  are  denoted  Pp^,  Cqs,  and  Aqj,  respectively.  Thus  Eq  (C.4) 
may  be  written 
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Equation  (C.8)  does  not  have  quite  the  same  form  as  Eq.  (B.2).  To 
get  the  desired  form,  we  Introduce  another  transformation,  this  time  of  the 
Independent  variable; 


Z-Z(s)  .  (CIO) 

By  the  chain  rule,  we  find  that  Eq.  (C.8)  can  be  written  as  follows; 
(dW/aZ)(dZ/ds)  -  I(p2VPosCo,5)/(p^2A^p^Co5)J  (p,/p^^Co55)WdW/at' .  (C.  1 1 ) 
Putting 

dz/as  =  V<hCo5^>''<Ps^ Vo'^o®)  fc. '  2) 


gives  Eq.  (C.  1 1 )  the  form  of  Eq.  (B.2): 


^  _  PsW  w 
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Equation  (C.  12)  must  now  be  solved  subject  to  the  condition  that  Z  =  0  when 
s  •  Sq.  The  required  solution  Is 


Z 


Sq  '  Ps  ^0  Po  ^0  ' 


(C.14) 


Z  Is  called  the  distortion  range  variable.  Equation  (C.  1 3)  has  the  desired 
form. 
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In  this  section  the  simplification  of  the  nonlinear  geometrical 
acoustics  equation,  Eg.  (B.10)  is  discussed.  The  approach  taken  to  reduce  the 
equation  is  similar  to  that  used  in  the  previous  chapter  to  simplify  the 
linear  geometrical  acoustics  equation,  Eq.  (3-C.17).  Since  self-refraction 
is  assumed  to  be  negligible,  the  eikonal  equation  from  linear  geometrical 
acoustics,  Eq.  (3-C.18),  is  used  to  simplify  the  nonlinear  geometrical 
acoustics  equation.  The  transport  equation  is  reduced  via  two 
transformations;  one  on  the  dependent  variable,  the  pressure,  and  the  other 
on  the  independent  variable,  the  path  length.  The  transformation  of  the 
pressure  corrects  for  the  geometrical  spreading  of  the  wavefront;  the 
transformation  of  the  path  length  corrects  for  the  increase,  or  decrease,  in 
distance  required  for  the  waveform  to  distort  a  prescribed  amount.  In  the 
finite  amplitude  case,  the  transforms  are  Eqs.  (C.9)  and  (C.  14).  The 
corresponding  transforms  in  the  linear  case  are  Eqs.  (3-E.7)  and  (3-E.  1 1 ), 
respectively.  Note  that  in  the  linear  case  the  transform  of  the  dependent 
variable  does  not  depend  on  the  static  density  or  sound  speed.  This  is 
because  in  the  linear  problem  the  level  of  approximation  is  only  first  order, 
whereas  in  this  chapter  the  level  of  approximation  is  consistently 
maintained  at  second  order 


Now  that  the  equation  for  propagation  in  an  inhomogeneous 
medium  has  been  cast  in  the  same  form  as  the  equation  for  plane  wave  in  a 
homogeneous  medium,  the  solution  may  readily  be  obtained.  The  plane  wave 
equation,  Eq.  (6.2),  is  satisfied  by  the  Eamshaw  solution  (see.  for  example. 


Blackstock  1972):  Given  the  boundary  condition 


P' =g(t)  =  g(t')  at  x  =  0  .  (D.l) 


the  solution  of  Eq.  (B.2)  Is 

II 

(D.2) 

Where 

♦  -  t'  ♦  pxg((l>)/CQ2  . 

(D.3) 

For  the  case  of  an  Inhomogeneous  medium  the  boundary  condition  is 

P'  =  g(t)  =  g(t')  at  s  =  .  (D.4) 

This  transforms  to 

W  =  g(t')  at  Z  =  0 
The  solution  is  therefore  the  following  equation 

W-g(^)  . 

where 

4>  =  t’  + 

E.  Simplification  of  the  Distortion  Range  Variable  Transformation 

In  this  section  the  integral  for  the  distortion  range  variable  Z, 

Eq.  (C.  1 4),  is  placed  in  the  more  compact  form  suggested  by  Morfey  ( 1 984a). 
Since  spreading  waves  In  a  mildly  Inhomogeneous  medium  are  similar  to 
spherical  spreading  waves,  Morfey’s  approach  Is  to  define  the  distortion 
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range  variable  for  an  Inhomogeneous  medium  in  a  form  similar  to  that  for  a 
spherically  spreading  wave  in  a  homogeneous  medium  (  Z  =  Sq  £n  (s/Sq)  ; 
Blackstock  1966).  liorfey's  definition  is  as  follows; 

Z  =  Sq  Jtn  (sG/Sp)  .  (E.l) 

where  s  is  the  path  length,  Sq  is  the  reference  path  length,  and  6  is  the 
dimensionless  distance  modification  factor.  As  is  shown  below,  6  embodies 
the  effects  of  both  ray  tube  geometry  and  ocean  inhomogeneity  upon  the 
finite  amplitude  distortion.  If  the  quantity  6  is  equal  to  1 ,  Eq.  (E.  1 )  reduces 
to  the  form  for  a  spherically  spreading  wave  in  a  homogeneous  ocean.  As 
Morfey  (1984a)  points  out,  if  G  is  greater  than  1,  finite  amplitude  distortion 
is  greater  than  it  would  be  in  a  homogeneous  medium.  The  converse  is  true 
if  6  is  less  than  1. 

Starting  with  Eq.  (Cl 2),  we  strive  to  force  it  into  the  form  of 
Eq.  (E.  1 ).  This  is  done  by  differentiating  both  equations  with  respect  to  the 
path  length  s,  and  then  equating  the  differentials.  First,  in  terms  of  a  new 
thermodynamic  variable  a, 

a  =  p(PoCo®)'''^  .  (E.2) 

Eq.  (C.  1 2)  may  be  expressed  as  follows: 

Z  -  (a/a,)(  V  *  .  (E-3) 

Where  a,  is  a  evaluated  at  the  source  position.  If  Eq.  (E.3)  is  differentiated 
with  respect  to  the  the  path  length  s,  we  are  left  with  the  integrand.  We 


define  a  new  variable  Y(s)  equal  to  the  Integrand: 

Y(s)  s  (a/a^XAo/Ao,)-’'^  .  (E.4) 

in  the  next  chapter,  we  use  a  computer  ray  model  based  on  the  assumption 
that  the  ocean  is  azimuthally  symmetric.  Using  this  assumption,  we 
substitute  for  the  ratio  of  the  ray  tube  areas,  Aq/Aq3,  using  a  relation  given 
by  Foreman  (1983), 

V^)9  “  <r^/So2)  (cos  6/  cos  6^)  ,  (E.5) 

where 

C  =  (dz/d0^),  ,  (E.6) 

r  is  the  radial  range,  z  is  the  depth,  and  6,  is  the  launch  angle.  By  defining  a 

new  variable  B(s),  we  see  that 

Y(s)  =  SoB(s)  ,  (E.7) 

where 

B(s)  s  (r^)"’'^(cos  6/  cos  6^)’’'^  (a/a^)  .  (E.8) 

Equation  (C.  12)  is  now  in  the  form  of  Eq.  (E.7). 

Now  differentiate  Eq.  (E.  1 )  with  respect  to  the  path  length  s  and 
equate  the  result  with  Eq.  (E.7).  The  differential  of  Eq.  (E.  1 )  is  as  follows: 

dZ/ds  =  (Sq/G)  (dG/ds)  ♦  Sq/s  .  (E.9) 

If  Eq.  (E.7)  is  used,  it  may  be  shown  that 


1/G  (dG/ds)  *B(s)-  1/s 


(E.IO) 
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Solving  Eq.  (E.  10)  for  G,  and  defining  a  new  variable  F(s),  we  see  that 


G  =  expF(s)  , 

(Eli) 

where 

F(s)  s  Jq  f(s)  ds  , 

(E.12) 

and 

f(s)sB(s)-l/s  . 

(E.13) 

Thus  the  distortion  range  variable  may  be  expressed  in  the  form 
of  Eq.  (E.  1 ),  and  the  quantity  G  may  be  evaluated  using  Eq.  (E.  1 1 )  if  both  the 
ray  path  and  the  environment  along  the  ray  path  are  known.  Thus  G 
embodies  both  the  effects  of  the  ray  tube  geometry  (via  the  area  In  Eq.  E.8) 
and  the  environment  (via  a  in  Eq.  £.8)  on  finite  amplitude  propagation. 

As  was  mentioned  in  Chapter  2,  the  main  difference  between 
linear  and  nonlinear  hydrodynamics  equations  can  be  explained  using  the 
mathematical  term  ranking  system.  To  form  the  linear  lossless 
hydrodynamics  equations  for  inhomogeneous  fluids,  we  retain  first-order 
terms  and  one  type  of  second-order  terms,  inhomogeneity  terms.  To  form 
the  nonlinear  lossless  hydrodynamics  equations,  we  include  one  more  set  of 
second-order  terms:  quadratic  nonlinear  terms.  These  terms  govern  the 
finite  amplitude  behavior. 
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CHAPTER  5 


NUMERICAL  EVALUATION  OF  WAVEFORMS 


A.  Introduction 

In  this  chapter  we  discuss  a  scheme  for  numerically  implementing 
weak  shock  theory  for  a  lossy  inhomogeneous  ocean.  Recall  that  in 
Chapter  4  the  transport  equation  for  finite  amplitude  signals  propagating 
through  an  inhomogeneous  ocean  was  reduced,  via  two  transforms,  to  the 
form  of  that  for  plane  waves  propagating  through  a  homogeneous  ocean. 
Pestorius  (1973)  developed  an  algorithm  for  numerically  implementing 
weak  shock  theory  in  the  case  of  plane  waves  propagating  through  a 
homogeneous  medium  in  a  pipe.  We  start  this  chapter  by  reviewing 
Pestorius's  algorithm.  We  then  describe  the  differences  between  his 
algorithm  and  the  algorithm  as  it  is  used  in  this  work;  related  details  of 
the  algorithm  are  examined.  Next  the  testing  of  the  algorithm  is 
discussed,  and  the  results  of  the  testing  are  presented.  In  connection  with 
the  testing,  typical  results  of  the  program  are  shown. 

Pestorius's  algorithm  involves  propagating  an  arbitrary  finite 
amplitude  waveform  a  small  distance  using  a  numerical  Implementation  of 
weak  shock  theory.  The  waveform  is  then  Fourier  transformed,  and  the 
viscosity  and  tube  wall  effects  incurred  over  that  small  distance  are 
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accounted  for  in  the  frequency  domain.  The  procedure  is  repeated  until  the 
desired  propagation  distance  has  been  reached.  Note  that  by  using 
Pestorius’s  algorithm  we  may  account  for  the  loss  terms  neglected  at  the 
beginning  of  the  analytic  development. 

Several  differences  exist  between  the  algorithm  as  set  forth  by 
Pestorius  and  the  modified  algorithm  used  in  this  study.  All  the 
differences  result  from  the  fact  that  Pestorius’s  algorithm  was  designed  to 
aid  in  the  study  of  nonlinear  propagation  of  plane  waves  through  a 
homogeneous  medium  in  a  pipe,  whereas  we  are  concerned  with  spreading 
waves  in  an  inhomogeneous  ocean.  Obviously  the  loss  mechanisms 
themselves  are  different;  we  account  for  the  viscosity  and  relaxation  of 
seawater.  From  a  computational  point  of  view,  two  major  differences  are 
( 1 )  that  our  absorption  coefficients  change  with  position,  and  (2)  that  we 
must  calculate  the  distance  over  which  the  absorption  acts.  Other 
differences  are  related  to  the  choice  of  step  size  and  starting  conditions. 
We  recalculate  the  step  size  depending  on  the  wave's  current  position, 
whereas  Pestorious  used  a  constant  step  size.  Pestorius  started  his 
waveform  at  range  zero;  we  must  select  a  reference  position. 

A  flowchart  of  the  modified  Pestorius  algorithm  is  shown  in 
Fig.  5.1.  The  modified  Pestorius  algorithm  is  implemented  in  the  program 
PLPROP 


FIGURE  5.1 

FLOWCHART  OF  THE  MODIFIED  PESTORIUS  ALGORITHM 


In  this  section  we  take  a  more  detailed  look  at  the  modified 
Pestorius  algorithm.  We  first  examine  the  weak  shock  subroutine 
WAVPROP,  and  then  focus  our  attention  on  the  calculation  of  attenuation  and 
dispersion.  The  section  closes  with  a  description  of  some  considerations  in 
the  numerical  evaluation  such  as  choice  of  nondimensional  variables,  step 
size,  and  starting  range. 

1.  The  Weak  Shock  Subroutine 
The  program  PLPROP  is  centered  around  the  weak  shock 
propagation  subroutine  WAVPROP.  This  subroutine  uses  the  Earnshaw 
solution  and  the  relative  shock  arrival  time  equation  (see,  for  example, 
Blackstock  1972)  infinite  difference  form  (Pestorius  1973,  Eqs.  4.2  and 
4.3).  From  the  Earnshaw  solution  which  describes  the  distortion  of  the 
continuous  portion  of  the  waveform,  one  finds  that  the  delay  time 
associated  with  a  given  wavelet  u  is 

t'Kk*  1  )hl  =  t'(kh)  -  pCo’^u  lt'(kh))h  ,  (B.  1 ) 

where  t'l(k+ 1  )h]  gives  the  position  of  the  wavelet  after  k+ 1  Incremental 
steps  of  size  h,  and  t'(kh)  gives  the  position  after  k  steps.  The  delay  time 
for  each  wavelet  is  calculated  individually  because  each  wavelet 
propagates  with  a  different  velocity.  The  finite  difference  form  of  the 
relative  shock  arrival  time  equation  gives  the  value  of  t'  associated  with  a 
particular  shock,  denoted  t,',  as  follows: 


Vl(k*  1  )hj  -  t;(kh)  -  phCo-2{u2lkh.t;(kh)l  ♦  u,lkh,t;(kh)l)/2 ,  (B.2) 
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where  t5'((k+ 1  )h]  is  the  value  of  tg’  after  k+ 1  steps  of  size  h,  and 
u,  zfkh.tg'fkh)]  is  the  value  of  particle  velocity  after  k  increments  of  size  h. 
These  two  equations.  Eqs.  (B.  1 )  and  (B.2),  are  all  that  is  required  to  describe 
the  propagation  of  weak  discontinuities  within  finite  amplitude  waves. 

Several  disadvantages  are  associated  with  having  to  transform 
back  and  forth  between  the  time  and  frequency  domains.  These  problems 
are  discussed  in  detail  by  Pestorius  (1973,  p.  91 ).  We  now  address  the  most 
significant  of  these  problems.  Equal  time  intervals  between  waveform 
points  are  required  by  the  fast  Fourier  transform  (FFT)  which  is  used  to 
transform  the  waveform  into  the  frequency  domain.  After  one  step  in  the 
the  WAVPROP  subroutine,  however,  the  waveform  points  are  no  longer 
spread  at  equal  time  intervals.  To  ensure  that  the  points  were  equally 
spread  through  time,  Pestorius  wrote  the  subroutine  RESAMP.  The 
subroutine  RESAMP  unshocks  the  waveform  by  spreading  the  shock  over  at 
least  one  time  interval.  Therein  lies  the  problem;  this  does  not  correspond 
with  the  physics  of  the  situation.  The  effects  of  the  problem  are  reduced 
by  having  many  data  points  closely  spaced  in  time.  As  can  be  seen  in  the 
later  section  on  testing,  the  effects  of  this  problem  are  small. 

Another  problem  is  aliasing  (see,  for  example,  Oppenheim  and 
Schafer  1975).  The  waveform's  spectrum  shows  some  degree  of  aliasing 
because  the  resampled  time  waveform  contains  very  high  frequency 
information.  The  spectra  of  the  waveforms  studied  in  this  report  have  a 
-6  dB/octave  slope,  and  the  effect  of  aliasing  is  therefore  small.  Also  the 


amount  of  attenuation  increases  with  frequency  and  acts  as  a  natural 
anti-aliasing  filter. 

2.  The  Application  of  Attenuation 

In  the  introduction  we  mentioned  that  differences  exist  between 
Pestorius's  algorithm  and  the  modified  algorithm  implemented  in  PLPROP. 

It  was  also  noted  that  the  major  differences  center  on  the  two  components 
required  to  calculate  the  attenuation:  the  absorption  coefficients  and  the 
distance  over  which  they  act.  Since  Pestorius’s  algorithm  was  intended  for 
plane  waves,  our  modified  version  operates  in  units  of  equivalent  plane 
wave  distance.  Consequently  the  true  distance  must  be  calculated  before 
the  attenuation  can  be  applied.  To  obtain  the  true  propagation  distance  we 
use  Eq.  (4-E.  1 ).  The  more  compact  form  of  the  distortion  range  variable 
requires  the  calculation  of  the  quantity  6. 

The  numerical  evaluation  of  G  requires  knowledge  of  both  the 
acoustic  ray  path  and  the  environment  along  the  ray  path.  Knowledge  of  the 
environment  is  usually  obtained  by  measurements  made  during  the  course  of 
an  experiment.  For  our  purposes  we  consulted  tabulated  data  and  chose 
plots  of  temperature  and  salinity  versus  depth  (hereafter  referred  to  as 
temperature  and  salinity  profiles)  typical  of  the  North  Atlantic  Ocean.  The 
profiles  are  shown  in  Figs.  5.2  and  5.3.  Once  the  temperature  and  salinity 
profiles  are  known,  the  sound  speed  profile  may  be  readily  calculated 
(Lovett  1978). 

As  was  pointed  out  in  Chapter  3,  the  acoustic  ray  paths  depend 
solely  on  the  sound  speed.  In  Chapter  4  it  was  noted  that  finite  amplitude 
waves  follow  the  same  ray  paths  as  their  small-signal  counterparts.  We 
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may  therefore  use  a  linear  acoustic  ray  model  to  calculate  the  ray  paths 
for  the  finite  amplitude  case;  the  ray  model  used  in  this  study  is  MEDUSA 
(Foreman  1983).  MEDUSA  assumes  an  azimuthally  symmetric  environment 
and  is  capable  of  calculating  the  ray  paths  even  if  the  sound  speed  varies 
with  both  depth  and  range.  However,  in  order  to  keep  our  analysis  simple, 
we  do  not  make  use  of  this  capability.  We  are  therefore  assuming  that  an 
ocean  may  be  modeled  by  a  single  sound  speed  profile,  a  stratified  ocean. 

The  RAYFAN  subroutine  in  htDUSA  was  modified  to  output 
information  about  the  ray  path.  The  information  includes  path  length, 
range,  depth,  cos  6,  and  i.  The  angle  6  is  defined  in  Fig.  3.2,  and  i  is  defined 
in  Eq.  (4-E.6).  The  cos  9  and  ^  values  are  required  to  calculate  the  ray  tube 
area  ratio  defined  in  Eq.  (4-E.5).  This  ray  path  information  is  calculated 
and  stored  by  RAYFAN  for  each  MEDUSA  step.  MEDUSA  varies  its  step  size  in 
order  to  maintain  a  constant  degree  of  accuracy.  It  is  therefore  assumed 
that  this  step  size  is  sufficiently  small  to  accurately  perform  the 
numerical  integration  outlined  below. 

Since  both  the  environmental  information  and  the  ray  path 
information  are  available,  the  quantity  G  may  be  numerically  evaluated 
Applying  a  second-order  numerical  integration  scheme  to  Eq.  (4-E.  1 1 ), 

Morfey  (1984a)  obtained  the  following  equation: 


f(s)  ds  »  h(f, 


♦f. 


\-\ 


)/2-hV(-fi'_,)/12 


(B.3) 


where  h  Is  the  step  size  (s,  -  s,.,)  in  arc  length.  The  function  f(s)  is  defined 


In  Eq.  (4-E.12),  and  f'(s)  is  the  derivative  of  f(s)  with  respect  to  s.  To 
obtain  the  quantity  G,  we  integrate  along  the  ray  path  using  Eq.  (6.3)  and 
then  substitute  the  resulting  value  into  Eq.  (4-E.  10).  The  program  CALCG 
implements  Eq.  (B.3)  and  calculates  the  quantity  G  at  each  of  the  MEDUSA 
steps.  Since  CALCG  requires  ray  trace  data  as  an  input,  MEDUSA  is  run 
first,  and  a  file  of  ray  trace  data  is  stored.  CALCG  appends  the  value  of  G 
for  each  step  to  the  ray  trace  data  file.  The  modified  ray  trace  data  file  is 
used  as  input  to  PLPROP.  In  this  way  PLPROP  has  access  not  only  to  the 
environmental  information  along  the  ray  path,  but  also  to  the  quantity  G. 
PLPROP  can  therefore  convert  the  distortion  range  variable  into  the  true 
distance  which  is  required  to  evaluate  the  amount  of  attenuation  associated 
with  a  propagation  step. 

The  second  major  difference  between  Pestorius’s  algorithm  and 
the  modified  algorithm  implemented  in  PLPROP  is  that  our  absorption 
coefficients  vary  along  the  propagation  path  whereas  Pestorius’s  were 
constant.  The  variation  of  the  temperature  and  salinity  causes  a 
corresponding  variation  in  the  absorption  coefficients. 

Empirical  relations  derived  by  Francois  and  Garrision  (1982)  are 
used  to  calculate  the  absorption  coefficients.  Francois  and  Garrision 
considered  three  sources  of  attenuation:  viscosity  and  the  relaxation 
mechanisms  associated  with  magnesium  sulfate  and  boric  acid.  Each  of  the 
relaxation  mechanisms  has  a  small  amount  of  dispersion  associated  with  it. 
The  dispersion  is  calculated  using  Blackstock's  technique  (1985).  The 
dispersion  due  to  each  of  the  two  relaxation  mechanisms  is  added. 

The  environmental  information  used  as  inputs  to  Francois  and 
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Garrision’s  empirical  relations  was  obtained  from  the  temperature  and 
salinity  profiles  via  a  cubic  T-spline  (see,  for  example,  Foreman  1983).  For 
the  splines  to  be  used,  the  depth  at  the  position  of  interest  must  be  known. 
The  depth  is,  however,  only  known  at  specific  points  along  the  ray  path 
since  MEDUSA  steps  through  range.  To  evaluate  the  depth  for  points  in  the 
middle  of  a  MEDUSA  step,  PLPROP  uses  a  linear  interpolation  scheme. 

Once  the  true  propagation  distance  and  the  absorption  coefficients 
have  been  calculated,  the  attenuation  and  dispersion  may  be  applied  to  the 
waveform.  However,  PLPROP  applies  the  attenuation  due  to  viscosity  only 
when  the  waveform  contains  no  shocks.  The  reason  for  this  is  that  weak 
shock  theory  implicitly  accounts  for  most  of  the  viscous  attenuation  when 
shocks  are  present.  Thus  if  PLPROP  were  to  apply  viscous  attenuation  to  a 
waveform  containing  shocks,  the  viscous  attenuation  would  have  been 
accounted  for  twice,  once  by  the  weak  shock  theory  and  once  by  PLPROP.  To 
separate  the  two  cases,  PLPROP  measures  the  distance  the  wave  propagates 
both  with  and  without  shocks,  and  then  applies  absorption  accordingly. 

The  reason  why  the  losses  associated  with  weak  shock  theory  are 
assumed  to  be  due  to  viscosity  is  as  follows:  In  weak  shock  theory  the 
losses  are  assumed  to  be  concentrated  at  the  shock.  Since  the  shock 
contains  mostly  high  frequency  information,  and  the  attenuation  mechanism 
which  is  dominant  at  high  frequencies  is  viscosity,  viscosity  alone  must  be 
responsible  for  the  losses  at  the  shock.  Other  attenuation  mechanisms, 
such  as  relaxation,  do  not  attenuate  the  shock  at  a  rate  sufficient  to  stop 
the  shock  from  becoming  multivalued. 


PLPROP  works  with  nondimensionalized  time  and  distance 
variables.  The  time  is  nondimensionalized  by  dividing  it  by  a  time  t^. 
characteristic  of  the  input  waveform.  For  example,  in  the  case  of  a  weak 
shock  with  an  exponentially  decaying  tail,  a  convenient  characteristic  time 
is  the  initial  1/e  decay  time.  The  characteristic  distance  used  to 
nondimensionalize  the  distortion  range  variable  is  defined  as  follows. 

R,  =  Cgt^/pe  .  (B.4) 

The  definition  for  R,  presented  in  Eq.  (B.4)  is  preferred  over  c^t^  alone 
because  Eq.  (B.4)  represents  a  distance  over  which  the  waveform  would 
undergo  a  significant  amount  of  distortion.  In  the  case  of  a  sinusoidal 
wave,  R,  is  the  shock  formation  distance.  PLPROP  propagates  the  waveform 
in  steps  of  nondimensional  distance  sigma  o,  defined  as  follows: 

0  =  z/  R.  .  (B.5) 

Several  factors  must  be  considered  in  choosing  PLPROP  s  step 
size,  denoted  Ad.  Pestorius  (1973)  tried  a  variety  of  step  sizes  and  settled 
on  the  value  of  0.01.  Pestorius  also  noted  that  several,  typically  10,  of  the 
0.01  Ad  steps  should  be  taken  before  the  attenuation  and  dispersion  are 
applied.  This  reduces  the  number  of  FFTs,  and  hence,  the  effect  of  the 
errors  associated  with  the  transformations  between  the  time  and 
frequency  domain. 

Step  size  considerations  for  the  modified  Pestorius  algorithm 
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implemented  in  PLPROP  are  somewhat  different  than  those  used  in  the 
original  algorithm.  One  of  the  consequences  of  Eq.  (4-E.  1 )  is  that  the  true 
distance  corresponding  to  0.01  Ac  increases  with  the  distance  the  wave  has 
propagated.  Close  to  the  source,  where  the  signal  is  stronger  and  finite 
amplitude  effects  are  important,  ten  0.01  Ac  steps  may  correspond  to  a 
very  small  distance.  To  calculate  the  number  of  0.01  Ac  steps  to  take, 
PLPROP  examines  the  absorption  coefficient,  calculated  in  dB/m,  from  the 
center  frequency  cell  of  the  FFT.  The  magnitude  of  the  reciprocal  of  this 
coefficent  gives,  for  that  particular  frequency,  the  propagation  distance 
required  for  a  1  dB  drop  in  the  amplitude.  The  waveform  is  then  propagated 
the  number  of  0.01  Ac  steps  corresponding  to  the  1  dB  drop  distance.  This 
differs  from  the  fixed  value  of  ten  0.01  Ac  steps  used  by  Pestorius. 
Compared  with  Pestorius's  original  approach,  the  number  of  applications  of 
absorption  near  the  beginning  of  the  ray  path  is  reduced.  This  corresponds 
to  the  decreased  role  of  absorption. 

On  the  other  hand,  the  number  of  0.01  Ac  steps  is  never  permitted 
to  go  below  10.  When  the  wave  is  far  from  the  source,  the  1  dB  drop 
distance  may  correspond  to  less  than  one  0.01  Ac  step.  Thus  to  combat  the 
errors  associated  with  transforming  between  the  time  and  frequency 
domains,  the  number  of  0.01  Ac  steps  between  applications  of  absorption  is 
fixed  at  10. 

Another  difference  between  the  algorithm  as  developed  by 
Pestorius  and  the  algorithm  used  in  PLPROP  is  the  requirement  for  a 
starting,  or  reference,  range.  Because  Pestorius  dealt  solely  with  plane 
waves,  his  starting  range  was  always  zero.  The  reference  range  used  in 


PLPROP  was  developed  by  Morfey  (1984)),  who  chose  the  reference  range  to 
be  equal  to  the  characteristic  distance  R,.  Morfey  formed  the  ratio  of  the 
nondimenslonal  shock  pressure,  AP/PqCq^,  to  the  nondimenslonal  time 
constant,  CqVR.  He  observed  that  the  ratio  Is  approximately  equal  to  1/p  if 
the  nondimenslonal  shock  pressure  is  0.06  and  the  range  Is  R«.  The  value  of 
0.06  for  the  ratio  AP/PqCq^  corresponds  to  the  upper  limit  on  weak  shock 
theory  (Pestorlus  and  Williams  1974).  Consequentially,  Morfey  chose  the 
model's  starting  conditions  to  correspond  to  a  nondimenslonal  pressure  of 
0.06  and  a  starting  range  of  R.. 

Another  feature  of  considerable  importance  is  PLPROP’s  ability  to 
selectively  Include  the  physical  mechanisms  that  affect  the  propagation. 

The  finite  amplitude  effects  may  be  switched  off  altogether  or  remain  in 
effect  only  until  the  wave  reaches  a  certain  range.  Viscosity  and 
relaxation  may  be  switched  on  Individually.  The  dispersion  due  to 
relaxation  may  be  Included  or  not.  The  absorption  coefficients  may  be 
calculated  from  the  local  environment  or  may  be  based  on  some  average 
value.  The  ability  to  switch  the  mechanisms  on  or  off  permits  each  of  them 
to  be  examined  individually. 

C.  Verification  of  the  Algorithm 

In  this  section  we  present  some  of  the  results  of  the  tests 
performed  to  verify  the  accuracy  of  the  modified  Pestorlus  algorithm.  The 
first  program  discussed  is  CALC6  which  calculates  the  quantity  G.  The 
second  program  discussed  Is  PLPROP.  In  each  case  the  numerical  results  of 
the  programs  are  compared  to  results  obtained  from  analytic  solutions. 
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The  percent  error  in  the  numerical  results  are  calculated. 

As  mentioned  earlier  in  this  chapter,  the  quantity  G  is  calculated 
using  Eq.  (B.3).  The  program  CALC6  evaluates  G  from  ray  path  information 
provided  by  the  ray  tracing  program  MEDUSA.  To  verify  the  accuracy  of 
CALCG,  we  compare  its  results  to  the  known  analytic  expression  for  G  for 
the  case  of  a  constant  gradient  sound  speed  profile  (Morfey  1984a).  Morfey 
has  shown  that,  in  this  case,  G  is  related  to  6  by  the  following  equation: 


sin  e  -  sin  AqX  /  tan  (n/4  ♦  e/2)  \ 


/  sin  0  -  sin  0o\  / 


0q)  cose/  \tan  (tt/4  +  0^/2) 


I 


(D.1) 


where  0  is  the  angle  between  the  ray  path  and  the  horizontal,  and  0q  is  the 
angle  0  at  the  reference  path  length.  In  our  test  the  dependence  of  sound 
speed  on  depth  was  as  follows; 


Co=  1462  *  0.01762  z 


(D2) 


where  z  is  the  depth  in  meters,  and  Cq  is  the  sound  speed  in  m/s.  The  sound 
speed  was  assumed  to  be  constant  with  range.  The  source  was  assumed  to 
be  at  a  depth  of  2500  m,  and  the  launch  angle  was  22.5®.  The  results  of 
Eq.  (D.l),  the  numerical  result  of  CALCG,  and  the  percent  error  are 
presented  in  Table  5.1.  As  can  be  seen  from  the  table,  the  error  is 
consistently  below  0.01  %. 

In  an  effort  to  check  the  accuracy  of  PLPROP,  we  compared  the 
program’s  results  with  results  obtained  from  the  analytical  expressions 
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developed  In  Appendix  C.  The  expressions  In  the  appendix  are  for  a  weak 
shock  with  an  exponentially  decaying  tall  which  Is  spreading  spherically 
through  a  homogeneous  medium.  The  expressions  give  the 
nondimenslonallzed  peak  pressure,  the  nondimenslonal  1/e  decay  time,  and 
the  nondimenslonal  relative  shock  arrival  time.  The  peak  pressure  Is 
normalized  to  the  Initial  peak  pressure.  The  1/e  decay  time  is  normalized 
to  the  Initial  1/e  decay  time.  The  shock  arrival  time  is  measured  relative 
to  the  travel  time  calculated  assuming  the  small-signal  sound  speed  Cq.  To 
make  the  comparison,  we  simulate  a  homogeneous  ocean  by  giving  MEDUSA  a 
constant  sound  speed  profile;  the  resulting  acoustic  ray  paths  are  straight 
lines.  If  one  of  these  ray  paths  Is  used  as  input  to  PLPROP,  and  the 
attenuation  In  PLPROP  Is  "turned  off",  PLPROP’s  results  should  be  the  same 
as  those  for  a  spherically  spreading  wave  In  a  lossless  homogeneous  ocean. 
A  comparison  of  the  analytic  and  numerical  results  Is  presented  In 
Table  5.2.  It  Is  seen  that  the  error  Is  generally  In  the  range  ±0.25  %. 

PLPROP  plots  the  time  waveforms  and  their  corresponding  energy 
spectra.  The  time  waveforms  and  energy  spectra  corresponding  to  the 
values  of  the  "numerical"  column  In  Table  5.2  are  shown  in  Figs.  5.4  and  5.5, 
respectively.  The  horizontal  axis  In  Fig.  5.4  is  in  units  of  nondimensional 
retarded  time.  The  vertical  axis  is  the  nondimensional  pressure.  The 
pressure  shown  here  and  that  listed  in  Table  5.2  is  the  equivalent  plane 
wave  pressure,  that  is,  the  purely  geometrical  effect  of  spherical 
spreading  has  been  removed.  The  plots  of  the  time  waveforms  from 
different  positions  along  the  ray  path  have  been  superimposed. 
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The  energy  spectra  in  Fig.  5.5  are  also  superimposed.  The 
horizontal  and  vertical  axis  in  Fig.  5.5  are  frequency  in  Hz  and  energy 
spectrum  level  in  dB// 1  J/m^-Hz.  The  energy  spectrum  levels  correspond 
to  the  energy  in  the  equivalent  plane  wave;  that  is,  the  effects  of  spherical 
spreading  are  removed.  The  frequency  scale  is  correct  as  shown.  The 
different  spectra  correspond  to  the  waveforms  from  different  points  along 
the  ray  path. 

The  above  results  indicate  that  the  weak  shock  propagation  in  the 
modified  Pestorius  algorithm  works  to  a  high  degree  of  accuracy,  at  least 
in  the  case  when  the  inhomogeneous  medium  is  made  to  appear  homogeneous. 
The  other  main  component  in  the  algorithm  is  the  application  of  absorption. 
It  was  tested  by  ( 1 )  checking  as  to  whether  the  absorption  coefficients  are 
correct,  and  (2)  propagating  a  sawtooth  wave  with  finite  amplitude  effects 
switched  off  and  the  absorption  switched  on.  The  resulting  waveform  was 
compared  with  that  obtained  from  a  separate  calculation  based  on  the 
Fourier  series  solution  of  the  sawtooth  wave.  Each  of  the  Fourier 
coefficients  was  attenuated  and  phase  shifted  to  simulate  the  application 
of  attenuation  and  dispersion  by  PLPROP.  The  results  of  both  calculations, 
the  Fourier  series  summation  and  PLPROP,  were  in  very  close  agreement. 

No  straightforward  way  of  verifying  the  accuracy  of  the 
combination  of  finite  amplitude  effects  and  absorption  exists.  However, 
since  the  two  major  components  of  the  program  work  well  separately,  it  is 
reasonable  to  assume  that  they  will  woric  well  together.  The  computer 
program  PLPROP  Is  presented  in  Appendix  D. 


CHAPTER  6 


RESULTS 


A.  Introduction 

In  this  chapter  the  effects  of  nonlinear  distortion  and  ordinary 
attenuation  and  dispersion  on  the  propagation  of  signals  In  an 
Inhomogeneous  ocean  are  Investigated.  The  investigation  is  conducted  using 
the  numerical  algorithm  discussed  In  the  previous  chapter.  Some  of  the 
results  have  beeh  reported  previously  (Cotaras,  Morfey,  and  Blackstock 
1984).  The  signals  examined  are  transients,  specifically  a  weak  shock  with 
an  exponentially  decaying  tall  (hereafter  referred  to  as  an  exponential 
pulse)  and  a  more  realistic  explosion  waveform  containing  one  bubble  pulse. 
First  the  ocean  environment  is  discussed,  and  two  specific  ray  paths  are 
selected.  Then  the  specific  waveforms  used  are  discussed  and  presented 
with  their  energy  spectra.  Next  the  effect  of  inhomogenelty  on  nonlinear 
distortion  Is  examined.  The  effect  Is  Investigated  by  propagating  an 
exponential  pulse  through  a  lossless  stratified  ocean  and  comparing  the 
result  with  that  obtained  for  a  lossless  homogeneous  ocean.  The  combined 
and  Individual  effects  of  nonlinear  propagation  and  ordinary  attenuation  and 


dispersion  are  then  examined.  These  effects  are  investigated  by 
propagating  the  same  exponential  pulse  through  a  stratified  ocean  and 
comparing  the  results  obtained  considering  ordinary  attenuation  and 
dispersion  only,  finite  amplitude  effects  only,  and  the  combination  of  the 
two.  We  also  briefly  examine,  in  terms  of  the  arrival  time  of  the  peak 
pressure,  the  role  of  dispersion  in  long  range  propagation.  Lastly  we  try  to 
answer  the  question,  "to  what  distance  are  finite  amplitude  effects 
Important?"  This  Is  done  by  comparing  and  contrasting  the  energy 
spectrum  obtained  considering  finite  amplitude  effects  over  the  entire 
propagation  path  with  those  obtained  by  neglecting  finite  amplitude  effects 
after  certain  distances,  namely  150  m  and  1 100  m. 

6.  Design  of  the  Numerical  Experiment 

1.  Bau^aths 

The  ocean  environment  and  the  ray  paths  selected  are  shown  in 
Fig.  6. 1 .  The  horizontal  axis  is  the  range  In  km,  and  the  vertical  axis  is  the 
depth  in  m.  As  mentioned  in  Chapter  5  the  ocean  is  assumed  to  be 
stratified.  The  salinity  and  temperature  profiles  used  to  calculate  the 
sound  speed  profile,  shown  in  Fig.  6.1,  are  the  same  ones  employed  In 
Chapter  5;  see  Figs.  5.2  and  5.3.  The  sound  speed  profile  is  shown  at  both 
0  km  and  75  km,  thereby  indicating  that  the  sound  speed  does  not  vary  with 
range.  The  two  ray  paths  shown  start  at  different  depths,  300  m  and 
4300  m,  but  have  the  same  launch  angle,  8®  down  from  the  horizontal.  (The 


FIGURE  6. 1 

OCEAN  ENVIRONMENT  AND  ACOUSTIC  RAY  PATHS 


angle  appears  to  be  greater  than  S'*  because  the  horizontal  and  vertical  axes 
have  different  scales.) 

The  two  particular  source  depths  were  selected  for  the  following 
reasons.  The  shallow  source  depth  was  chosen  because  ocean  acoustic 
measurements  are  commonly  made  at  this  depth.  The  deep  source  depth  was 
chosen  primarily  to  permit  the  study  of  an  explosion  waveform  that 
Includes  the  first  bubble  pulse.  More  details  are  given  in  the  next  section. 

The  criteria  for  selecting  the  ray  path  from  each  of  the  two 
depths  are  the  same.  Because  reflections  and  caustics  cannot  be 
accommodated  by  our  computer  algorithm,  they  must  be  avoided.  At  the 
same  time,  a  long  ray  path  is  needed  in  order  to  determine  whether 
nonlinear  distortion  continues  to  accumulate  after  long  propagation 
distances.  The  rays  selected  travel  a  moderate  distance.  50  to  70  km. 
without  interacting  with  the  surface  or  the  bottom  and  without  passing 
through  a  caustic.  Shown  In  Figs.  6.2  and  6.3  are  the  families  of  rays  from 
which  the  rays  shown  in  Fig.  6.  l  were  picked.  In  Fig.  6.2  (the  shallow 
source)  the  ray  launch  angles  vary  in  0.5®  increments  from  0®  to  15® 
measured  down  from  the  horizontal.  The  ray  selected,  noted  in  the  figure, 
avoids  the  caustic  region  that  starts  at  about  33  km  on  the  upper  edge  of 
the  family  of  rays.  The  deep  source  family  (Fig.  6.3)  has  ray  launch  angles 
that  vary  in  0.5®  increments  from  -5®  to  10®.  The  ray  paths  all  reflect  from 
the  surface  before  passing  through  a  caustic  at  about  62  km.  The  ray 
selected  from  this  family  is  noted.  From  now  on  the  ray  path  starting  at 
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FIGURE  6.3 

FAMILY  OF  RAYS  FROM  5*  TO  10*  AT  0.5*  INCREMENTS  STARTING  FROM  4300  m 
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the  shallow  source  is  referred  to,  for  short,  as  the  shallow  path;  similarly 


the  other  path  is  referred  to  as  the  deep  path. 


2.  Signals 

The  basic  criteria  for  selecting  the  waveforms  to  study  are  as 
follows:  ( 1 )  They  must  be  waveforms  for  which  finite  amplitude  effects 
are  important.  (2)  They  must  be  as  realistic  as  possible,  and  yet  short 
enough  to  be  properly  sampled  by  our  computer  program.  The  signals 
examined  are  transient  pulses  similar  to  those  caused  by  underwater 
explosions.  The  explosives  commonly  used  in  long  range  propagation  studies 
generate  the  high  sound  pressure  levels  at  which  nonlinear  effects  play  an 
Important  role  in  propagation,  it  is  rare  that  a  long  range  propagation 
experiment  is  conducted  using  an  intense  periodic  source.  For  this  reason 
we  exclude  the  study  of  periodic  signals. 

With  regard  to  realism  of  the  waveforms  used,  It  is  noted  that 
Wakeley  (1977)  developed  an  empirical  relation  for  an  underwater  explosion 
waveform  which  includes  the  first  four  bubble  pulses.  The  waveform  has 
been  shown  to  give  a  close  fit  to  experimental  data;  it  is,  however, 
somewhat  complicated.  Morfey  (1985)  developed  a  simplified  version  of 
Wakeley  s  waveform  which  is  shorter  in  time  duration  and  is  therefore 
more  suitable  for  use  in  our  computer  program.  (The  simplified  Wakeley 
waveform  Is  shown  on  p.  III.)  The  approach  taken  by  Morfey  was  to  ( i ) 
truncate  the  Wakeley  waveform  at  the  zero  crossing  after  the  first  bubble 
pulse  and  (2)  remove  the  explicit  dependence  on  the  charge  weight  by 
selecting  a  nondimensional  time  base  that  is  related  to  the  charge  weight. 
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As  the  characteristic  time  used  for  nondimensionalization,  Morfey  picked 
the  initial  1/e  decay  time.  He  found  that  (Morfey  1984b)  use  of  the 
empirical  scaling  law  for  the  peak  pressure  (Arons  1954)  and  the  starting 
conditions  mentioned  in  Chapter  5  (nondimensional  shock  pressure  of  0.06 
at  a  range  of  R*)  leads  to  the  following  expression  for  the  1/e  decay  time; 

0^  =  59W’'3  ,  (6.1) 

where  0^  is  the  1/e  decay  time  in  ps  and  W  is  the  equivalent  TNT  charge 
weight  in  kg.  Equation  (6. 1 )  connects  the  initial  1  /e  decay  time  of  the  pulse 
with  the  charge  weight.  By  simplifying  Wakeley's  expression  so  that  it 
includes  only  one  bubble  pulse  and  by  using  0^  to  calculate  the 
nondimensional  time  T,  Morfey  (1985)  obtained  the  following  expression: 

P  -  exp(-T)  +  0.16exp[-(T  -  Tg)/T,]  -  (TT/2Tg)(l  +0.16T,)lsin(nT/Tg)|  ,  (6.2) 

where  Tg  and  T ^  are  the  nondimensional  bubble  pulse  period  and  bubble  pulse 
time  constant,  respectively.  The  bubble  pulse  period  is  the  time  between 
the  initial  pulse  and  the  bubble  pulse.  The  bubble  pulse  time  constant  is 
related  to  the  rate  of  rise  and  decay  of  the  bubble  pulse.  The  bubble  pulse 
period  and  time  constant  are  defined  in  the  following  expressions: 


Tg  =  35600  (z+  10.1)-5'6  , 

(6.4) 

T,  =  165(zM0.1)-’'2  , 

(6.5) 

where  z  is  the  depth  in  m. 


It  can  be  seen  from  Eq.  (6.4)  that,  for  large  depths,  the  bubble 
pulse  period  of  an  explosion  is  approximately  proportional  to  the  inverse  of 
the  depth.  Accordingly  the  explosion  waveform  from  a  source  at  300  m  has 
a  long  bubble  pulse  period,  almost  300  times  the  initial  l  /e  decay  time,  in 
order  to  maintain  the  accuracy  stated  in  Chapter  5,  our  computer  program 
requires  at  least  32  points  in  the  initial  l/e  decay  time.  Thus  9600  points 
are  required  if  the  shallow  source  explosion  waveform  is  to  include  the 
first  bubble  pulse.  The  total  number  of  points  in  our  time  waveform  is, 
however,  restricted  to  4096  by  computer  space  limitations.  Hence  we 
cannot  include  the  bubble  pulse  in  the  explosion  waveform  from  the  shallow 
source.  A  much  deeper  source,  however,  has  the  shorter  bubble  pulse 
period  which  fits  within  our  limitations.  We  chose  our  second  source  depth 
to  be  4300  m. 

A  simpler  waveform  is  more  appropriate  for  the  shallow  path 
since  the  details  of  the  more  complicated  waveform  cannot  be  included. 
Shown  in  Fig.  6.4  are  the  time  waveforms  and  corresponding  energy  spectra 
of  the  modified  Wakeley  waveform  for  the  300  m  source  (Fig.  6.4a)  and  the 
exponential  pulse  (Fig.  6.4b).  Both  waveforms  correspond  to  a  0.81 8  kg  TNT 
explosion  at  the  reference  range  of  0.4  m.  The  reference  range  is 
calculated  using  Eq.  (5-B.4)  which  assumes  a  peak  pressure  level  of  282.6 
dB// 1  pPa.  The  initial  l/e  decay  time  is  55  ps.  Since  the  time  waveforms 
and  frequency  spectra  of  Figs.  6.4(a)  and  6.4(b)  are  very  similar,  it  was 
decided  that  the  simpler  signal,  the  exponential  pulse  Fig.  6.4(b),  would  be 
the  signal  used  for  the  shallow  path. 


Shown  In  Figs.  6.5(a)  and  6.5(b)  are  the  modified  Wakeley 
waveforms  and  corresponding  energy  spectra  for  the  0.818  kg  and  22  7  kg 
TNT  explosions  at  a  depth  of  4300  m.  The  waveforms  are  calculated  at  the 
reference  ranges  of  0.4  m  and  1.1  m,  respectively.  Because  the  explicit 
dependence  on  charge  weight  has  been  removed  from  the  expression  for  the 
explosion  waveform,  Eq.  (6.2),  the  nondimenslonal  time  waveforms  are 
Identical,  if  the  time  waveforms  had  been  plotted  in  terms  of  dimensional 
units,  as  are  the  frequency  spectra,  the  waveforms  would  appear  different. 
The  larger  explosion  would  have  a  longer  bubble  pulse  period.  This  can  be 
seen  by  noting  that  the  22.7  kg  TNT  explosion  energy  spectrum  Is  shifted 
down  In  frequency  in  comparison  to  the  other  spectrum.  The  frequency 
shift  Is  caused  by  the  Increase  In  charge  weight.  Since  ordinary  absorption 
and  finite  amplitude  effects  scale  differently  with  frequency.  It  Is 
expected  that,  after  the  signals  propagate  a  moderate  distance,  their 
spectral  shapes  will  be  different.  It  is  for  this  reason  that  two  different 
charge  weights  are  considered. 


C.  Effect  of  Medium  Inhomogeneity  on  Nonlinearity 

The  effect  of  ocean  inhomogenelty  upon  finite  amplitude 
propagation  Is  small.  This  conclusion  was  first  quantitatively 
demonstrated  by  Morfey  (1984a).  To  understand  why  the  effect  is  small, 
recall  the  definition  of  the  distortion  range  variable,  Eq.  (4-E.  I ).  All  the 
effects  of  ocean  Inhomogenelty  are  embodied  In  the  nondimenslonal 
parameter  6.  The  value  of  (?  for  the  shallow  path  is  approximately  3. 
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Physically  this  means  that  to  undergo  the  same  amount  of  distortion  In  a 
homogeneous  ocean,  the  signal  would  have  to  propagate  three  times  as  far. 
This  is  not,  however,  a  large  enhancement  because,  as  can  be  seen  from 
Eq.  (4-E.  1 ),  the  distortion  range  varlable—and  hence  the  finite  amplitude 
effects— are  a  function  of  the  logarithm  of  the  distance  propagated. 

As  an  example  of  how  small  an  Influence  Inhomogenelty  has. 
consider  a  signal  propagating  along  the  shallow  path.  If  the  reference  path 
length  Is  1  m,  the  path  length  Is  20  km,  and  the  value  of  6  Is  3,  the 
distortion  range  variable  is  1 1.  If  the  ocean  is  homogeneous,  the  distortion 
range  variable  for  a  20  km  path  and  a  l  m  reference  range  Is  9.9.  Since  the 
amount  of  finite  amplitude  distortion  Is  proportional  to  the  distortion 
range  variable,  the  inhomogenelty  enhances  the  distortion  by  only  1  \%  over 
the  20  km  path  length. 

The  enhancement  of  finite  amplitude  effects  along  the  deep  path  Is 
even  smaller.  There  the  maximum  value  of  G  Is  approximately  l.  Hence  In 
terms  of  nonlinear  distortion  there  Is  effectively  no  difference  between 
propagating  through  a  homogeneous  ocean  and  propagating  along  the  deep 
path  through  an  Inhomogeneous  ocean.  It  Is  therefore  concluded  that  the 
effect  of  inhomogenelty  upon  nonlinear  distortion  is  small. 

We  now  quantitatively  demonstrate  the  small  effect  of 
Inhomogenelty  on  distortion  by  comparing  the  waveform  for  a  homogeneous 
ocean  with  that  for  an  inhomogeneous  ocean;  examine  Fig.  6.6.  The  initial 
waveform  is  the  exponential  pulse  shown  in  Fig.  6.4(b)  but  on  an  expanded 
time  scale.  The  other  waveforms  result  from  ( 1 )  propagating  the  signal 
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58. 1  km  along  the  shallow  path  through  a  lossless  ocean  and  (2)  propagating 
the  same  signal  to  the  same  distance  through  a  lossless  homogeneous  ocean. 
Note  that  the  strong  effect  of  geometrical  spreading  has  been  suppressed 
by  plotting  the  transformed  pressure  given  In  Eq  (4-C.9) 
nondlmenslonallzed  by  its  initial  value.  The  differences  In  amplitude  and 
relative  shock  arrival  time  between  the  two  resultant  waveforms  are 
clearly  small  (recall  that  Jin  G  is  small  compared  to  Jin  (s/Sq)). 

On  the  other  hand  nonlinear  effects  are  still  important;  reexamine 
Fig.  6.6.  The  shock  arrival  time  is  approximately  165  ps  (three  times  the 
Initial  1/e  decay  time)  In  advance  of  the  linear  theory  prediction,  and  the 
peak  pressure  is  approximately  1/3  of  Its  original  value.  The  nonlinear 
distortion  may  therefore  be  thought  of  as  "stretching”  the  initial  signal  as 
well  as  attenuating  it.  Note  that  the  "lossless  ocean"  is  not  really  lossless 
because  the  nonlinear  distortion  takes  account  of  losses  at  the  shock. 
However,  losses  due  to  the  medium  which  affect  both  continuous  and 
discontinuous  waves,  such  as  relaxation,  are  not  explicitly  accounted  for; 
hence  the  term  lossless  ocean. 

0.  The  combined  Effects  of  Nonlinearity  and  Absorption 

Throughout  the  propagation  of  a  finite  amplitude  wave,  the  effects 
of  both  nonlinearity  and  ordinary  attenuation  and  dispersion  are  at  play.  To 
more  clearly  examine  the  role  of  each  in  the  propagation  process,  we  first 
examine  them  separately,  and  then  together.  Figure  6.7  shows  four 
waveforms,  one  of  which  Is  the  initial  waveform,  the  exponential  pulse  of 
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Fig.  6.4(b)  on  an  expanded  time  scale.  The  others  are  the  resulting 
waveforms  after  propagating  58. 1  km  along  the  shallow  path  considering 
( 1 )  ordinary  attenuation  and  dispersion  only,  (2)  finite  amplitude  effects 
only,  and  (3)  both  finite  amplitude  effects  and  ordinary  attenuation  and 
dispersion.  It  is  clear  from  the  figure  that  by  itself  neither  finite 
amplitude  effects  nor  ordinary  absorption  can  correctly  account  for  the 
shape  of  the  resulting  waveform;  both  are  required.  The  finite  amplitude 
effects  attenuate  the  wave  while  causing  it  to  try  to  form,  or  to  maintain, 
a  shock.  Ordinary  absorption  attenuates  the  wave,  thus  causing  the  wave  to 
become  rounded. 

It  is  interesting  to  examine  the  effect  of  dispersion  on  the 

position  of  the  peak  pressure.  Consider  the  propagation,  neglecting  finite 

amplitude  effects,  of  the  exponential  pulse  along  the  shallow  path.  The 

waveforms  in  Figs.  6.8(a)  and  6.8(b)  are  obtained  at  various  positions  along 

the  shallow  path,  from  the  reference  range  0.4  m  out  to  58.  l  km.  The 

calculations  were  made  first  without  dispersion  (Fig.  6.8a),  and  then  with 

dispersion  (Fig.  6.8b).  At  longer  ranges  the  effect  of  dispersion  is  clear;  it 

shifts  the  waveform  forward.  In  Fig.  6.8(a)  the  peak  pressure  continuously 

moves  backward;  i.e.,  the  effective  propagation  speed  of  the  peak  is  less 

than  Cq  The  same  is  true  of  the  waveforms  in  Fig.  6.8(b)  for  distances  up  to 

8.6  km,  but  beyond  that  distance  dispersion  pulls  the  peak  forward.  By  the 

* 

time  the  wave  has  reached  the  maximum  distance  of  58. 1  km,  the  shift  of 
the  peak  pressure  due  to  dispersion  is  approximately  0.8  of  the  initial  l/e 
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decay  time.  This  would  be  a  significant  amount  of  time  if  one  were  trying 
to  add  signals  coherently. 

In  Fig.  6.9  we  again  examine  the  effect  of  dispersion,  but  this  time 
finite  amplitude  effects  are  also  Included.  In  this  case  the  waveforms  at 
the  lower  ranges  are  resolved.  Notice  that  in  Fig.  6.9(a)  the  peak  pressure 
moves  forward  until  a  range  of  3300  m  is  reached  and  then  moves 
backwards.  In  Fig.  6.9(b),  however,  the  forward  movement  of  the  the  peak 
pressure  is  monotonic.  The  forward  movement  is  due  first  to  finite 
amplitude  efferts  and  in  the  end  due  to  dispersion. 

Figure  6.9(a)  enables  us  to  examine  the  combined  effects  of 
nonlinearity  and  attenuation.  As  noted  above,  the  peak  pressure  in 
Fig.  6.9(a)  stops  moving  forward  as  the  signal  propagates  beyond  3300  m. 
One  might  interpret  this  to  Indicate  a  change  in  the  importance  of  finite 
amplitude  effects  relative  to  attenuation.  From  3300  m  on,  ordinary 
attenuation  is  the  dominant  mechanism  of  diminution,  whereas  for  ranges 
less  than  3300  m  finite  amplitude  effects  are  the  principal  mechanism.  It 
is  noted,  however,  that  the  amount  of  attenuation  beyond  3300  m  indicated 
in  Fig.  6.8(a)  is  less  than  that  over  the  corresponding  distance  in  Fig.  6.9(a). 
It  is  therefore  concluded  that,  even  though  finite  amplitude  effects  are  not 
dominant  beyond  3300  m,  they  are  still  noticeable.  Although  the  transition 
point,  3300  m,  applies  only  to  this  particular  example,  one  can  expect  a 
similar  behavior  for  other  signals  of  similar  initial  shape. 
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In  this  section  we  find  that  for  a  specific  frequency  range  and 
source  strength,  finite  amplitude  effects  may  be  neglected  beyond  a  certain 
distance.  Ordinary  absorption  appears  to  dominate  the  propagation  beyond 
that  distance.  One  explanation  for  the  importance  of  ordinary  absorption 
over  finite  amplitude  effects  is  that  the  effects  of  ordinary  absorption 
increase  as  exp(-as),  whereas  most  finite  amplitude  effects  depend  on  a 
(Eg.  5-B.5)  which  is  proportional  to  itn  (s/Sq). 

To  answer  quantitatively  the  question  posed  in  the  heading  of  this 
section,  we  conducted  a  numerical  experiment  involving  two  different 
explosion  pulses.  The  first  is  for  a  0.818  kg  TNT  explosioh  at  a  depth  of 
4300  m,  and  the  second,  a  22.7  kg  TNT  explosion  at  the  same  depth.  The 
source  waveforms  and  their  respective  spectra  are  shown  in  Figs.  6  5(a) 
and  6.5(b).  The  procedure  used  is  a  simple  one.  The  explosion  waves  are 
numerically  propagated  to  a  distance  of  23  km  along  the  deep  path 
accounting  for  attenuation  and  dispersion  over  the  entire  23  km  and 
accounting  for  nonlinear  effects  as  indicated  below. 

Case  A;  nonlinear  effects  neglected  entirely. 

Case  B:  nonlinear  effects  included  only  up  to  range  150  m, 

Case  C:  nonlinear  effects  included  only  up  to  range  1 100  m. 

Case  D:  nonlinear  effects  included  for  the  entire  23  km. 

Case  D  is  used  as  the  basis  of  comparison. 

Because  the  effective  duration  of  the  explosion  pulses  in  Fig.  6.5 
Is  much  greater  than  that  of  the  simple  exponential  pulses  (compare  the 


waveforms  of  Fig.  6.5  with  the  initial  waveform  of  p  ig  6  6),  time  waveform 
resolution  of  the  sort  shown  in  Figs  6  6  through  6  9  is  not  possiDie 
Interesting  results  may,  however,  be  foi^  by  comparing  the  spectra  of  the 
signals.  The  results  are  therefore  presented  in  the  form  of  energy 
spectrum  plots.  A  few  time  waveforms  are  shown  for  clarity  in 
Interpreting  changes  in  spectra 


The  results  of  the  numerical  experiment  involving  the  0  818  kg 
TNT  explosive  are  presented  in  Figs  6. 10  through  6  1 2  The  dotted  curve  in 
Fig.  6.10  is  the  Initial  energy  spectrum  at  the  reference  range,  0  4  m  The 
dashed  curve  Is  the  energy  spectrum  calculated  neglecting  finite  amplitude 
effects  over  the  entire  23  km.  Case  A.  The  solid  curve  in  the  figures  is  the 
energy  spectrum  calculated  accounting  for  finite  amplitude  effects  over 
the  entire  path.  Case  D.  Figure  6. 1 1  shows  a  comparison  of  Cases  B  and  D 
and  Fig.  6. 1 2,  Cases  c  and  D. 

We  start  our  discussion  by  examining  the  behavior  of  the  solid 
curve  (Case  D)  and  the  dashed  curve  (Case  A)  in  Fig.  6. 1 0.  Note  that,  as  for 
the  time  waveforms,  the  effect  of  geometric  spreading  has  been  removed. 
The  differences  between  the  two  curves  are  itemized  below. 

(a)  At  the  high  frequency  end  (above  15  kHz)  the  solid  curve  Is  higher. 

(b)  In  the  middle  range  (approximately  1.5-  15  kHz)  the  dashed  curve  is 
higher. 

(c)  At  the  low  frequency  end  (below  l  .5  kHz)  the  envelopes  of  the  two 
curves  are  about  the  same. 
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ENERGY  SPECTRA  OF  0.818  leg  TNT  EXPLOSION 
( I )  AT  THE  REFERENCE  RANGE  (0.4  m).  (2)  AFTER  PROPAGATING  23  km 
NEGLECTING  FINITE  AHPLITUDE  EFFECTS  (CASE  A).  AND  (3)  AFTER 
PROPAGATING  23  km  CONSIDERING  FINITE  AMPLITUDE  EFFECTS  (CASE  D) 


(d)  As  the  frequency  decreases  from  about  400  Hz.  the  solid  curve 
rises  above  the  dashed  curve. 

(e)  In  the  low  frequency  region  the  spectral  peaks  of  the  solid  curve 
occur  at  slightly  lower  frequencies  than  the  peaks  of  the  dashed 
curve. 

The  differences  cited  above  can  be  explained  in  terms  of  the 
nonlinear  distortion  of  the  wave.  Compare  the  two  time  waveforms 
Inserts.  Difference  (a)  is  probably  caused  by  the  steepening  of  the 
compression  part  of  the  bubble  pulse.  Difference  (b)  is  probably  due  to  the 
Increase  in  the  decay  time  of  the  first  peak  (as  the  shock  pulls  ahead  of  the 
first  zero).  Differences  (d)  and  (e)  are  probably  due  to  the  stretching  of  the 
time  Interval  between  the  initial  peak  and  the  bubble  pulse  peak. 

We  now  examine  the  Case  B  and  Case  C  spectra.  In  the  following 
discussion  of  Figs.  6. 1 1  and  6. 1 2  we  attempt  to  answer  the  question,  "to 
what  distance  are  finite  amplitude  effects  important?"  The  inclusion  of 
finite  amplitude  effects  up  to  150  m  (Fig.  6. 1 1 )  gives  a  23  km  spectrum  that 
follows  the  Case  D  spectrum  up  to  about  5  kHz.  Less  important  effects  are 
seen  at  the  very  low  frequencies.  Figure  6. 12  shows  that  the  Case  C 
spectrum  also  follows  the  Case  D  spectrum  up  to  about  6  kHz.  although  it 
does  so  more  closely  than  the  Case  B  spectrum.  The  tentative  conclusion, 
namely  that  nonlinear  effects  may  be  ignored  beyond  a  certain  distance, 
depends  on  the  frequency  range  in  which  one  is  Interested. 


We  now  turn  our  attention  to  the  case  of  a  22.7  kg  TNT  explosion. 
As  can  be  seen  In  Figs.  6.5(a)  and  6.5(b),  the  energy  spectra  of  the  two 
explosions  are  very  similar  except  that  the  larger  explosion  has  an  overall 
higher  spectrum  level  and  Is  shifted  down  in  frequency.  To  avoid  peak 
pressures  too  high  to  be  handled  correctly  by  weak  shock  theory,  we 
Increase  our  reference  range  to  1.1  m  (for  the  22.7  kg  pulse  only),  thereby 
making  the  peak  pressure  at  the  reference  range  the  same  for  both  the 
0.818  kg  and  the  22.7  kg  TNT  explosion  (282.6  dB// 1  pPa).  Because  of  the 
frequency  shift  of  the  spectrum,  the  nonlinear  distortion  of  the  two  pulses 
Is  much  the  same  (finite  amplitude  effects  scale  with  frequency).  The  only 
real  difference  Is  the  effect  of  attenuation,  which  should  be  less  for  the 
pulse  from  the  larger  charge. 

The  results  In  the  form  of  energy  spectra  are  presented  In 
Figs.  6. 1 3  through  6. 1 5.  The  Initial  spectrum  Is  shown  as  the  dotted  curve 
In  Fig.  6.13.  In  this  figure  the  dashed  curve  Is  the  energy  spectrum 
calculated  neglecting  finite  amplitude  effects  entirely.  Case  A.  The  solid 
curve  In  Figs.  6. 13  through  6. 15  Is  the  23  km  energy  spectrum  calculated 
accounting  for  finite  amplitude  effects  over  the  entire  path.  Case  D.  In 
Fig.  6.14  the  dashed  curve  Is  the  spectrum  for  Case  B;  In  Fig.  6. 15  the  dashed 
curve  Is  the  spectrum  for  Case  c. 

The  general  observations  made  of  the  solid  and  dashed  curves  In 
Fig.  6.13  are  the  same  as  those  made  of  Fig.  6. 10,  except  that  the 
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frequencies  cited  previously  are  higher.  The  differences  between  the  solid 
and  dashed  curves  of  Fig.  6.13  may  be  summarized  as  follows. 

(a)  Above  approximately  10  kHz  the  solid  curve  is  higher. 

(b)  From  about  500  Hz  to  10  kHz  the  dashed  line  is  higher. 

(c)  The  solid  curve  Is  slightly  higher  from  50  to  500  Hz. 

(d)  Over  the  same  frequency  range  as  (c)  the  spectral  peaks  of  the 
solid  curve  exhibit  a  slight  downward  shift. 

The  physical  explanations  for  these  differences  are  the  same  as  those  for 
the  differences  witnessed  in  Fig.  6. 10.  The  nonlinear  steepening  of  the 
compresslonal  part  of  the  bubble  pulse  accounts  for  difference  (a). 
Difference  (b)  Is  due  to  the  Increase  In  the  decay  time  of  the  first  peak. 
Differences  (c)  and  (d)  are  due  to  the  slight  increase  In  the  period  between 
the  Initial  pulse  and  the  bubble  pulse. 

We  now  address  the  question,  "to  what  distance  are  finite 
amplitude  effects  important?’  As  in  the  case  of  the  0.818  kg  explosion,  the 
answer  may  be  obtained  by  comparing  the  Case  8  and  Case  C  energy  spectra 
with  the  Case  D  spectrum.  The  Case  6  and  Case  C  spectra  (Figs.  6. 1 4  and 
6.15,  respectively)  follow  the  Case  D  spectrum  closely  for  frequencies 
below  4  kHz,  although  the  Case  C  spectrum  more  closely  duplicates  that  of 
the  Case  D.  It  Is  therefore  concluded  that  the  finite  amplitude  effects  can 
be  neglected  after  a  certain  distance,  and  that  the  distance  depends  on  both 
frequency  and  source  strength.  Since  the  differences  between  the  Case  C 
and  Case  D  spectra  are  approximately  the  same  for  both  the  0.818  kg  and 


22.7  kg  TNT  explosions,  It  is  also  concluded  that  the  differences  between 
the  two  spectra  due  to  absorption  are  small. 

In  summary  it  is  noted  that  the  prevailing  sentiment  that 
"nonlinear  effects  are  important  only  close  to  the  source"  Is  confirmed 
quantitatively.  Beyond  a  distance  that  depends  on  frequency  and  charge 
weight,  nonlinear  effects  may  be  neglected.  Notice,  however,  that  the 
calculations  are  for  rays  that  encounter  neither  reflections  nor  caustics 


CHAPTER  7 


SUMMARY 


In  this  report  the  propagation  of  finite  amplitude  signals  through 
an  inhomogeneous  ocean  is  investigated  both  analytically  and  numerically. 
The  theory  used  is  that  of  nonlinear  geometrical  acoustics.  The  amplitude 
of  the  signal  is  assumed  to  be  small  enough  that  self-refraction  may  be 
neglected.  Another  assumption  is  that  the  acoustic  field  consists  only  of 
outgoing  waves.  The  effects  of  reflections  and  focusing  are  not  considered. 
Losses  are  accounted  for  directly  in  the  numerical  routine;  the  analytic 
worK  is  performed  neglecting  all  losses  except  those  at  the  shock,  in  the 
numerical  study  the  ocean  is  assumed  to  be  stratified,  whereas  the  analytic 
work  is  for  a  fully  inhomogeneous  ocean  in  which  temperature,  salinity,  and 
density  vary  with  position. 

The  analysis  starts  with  the  presentation  of  a  scheme  for  ranking 
the  terms  in  the  hydrodynamics  equations  according  to  their  degree  of 
smallness.  The  ranking  scheme  is  then  used  to  simplify  the  hydrodynamics 
equations  for  lossless  inhomogeneous  fluids  for  (i)  small-signal  and 


(2)  finite  amplitude  waves.  From  these  simplified  equations  the  theories  of 
linear  and  nonlinear  geometrical  acoustics  are  developed.  The  eikonal 
equation  is  obtained  and  found  to  be  the  same  for  both  small-signal  and 
finite  amplitude  waves.  The  transport  equation  is,  however,  different  for 
the  two  cases.  An  equation  for  the  ray  paths  is  derived  from  the  eikonal 
equation  and  the  relations  of  vector  calculus,  in  the  case  of  a  constant 
gradient  sound  speed  profile,  the  ray  paths  are  found  to  be  circular  arcs. 
The  two  transport  equations  are  found  to  be  equivalent  to  the  first-order 
progressive  wave  equations  for  small-signal  and  finite  amplitude  waves, 
respectively.  All  the  analysis  is  carried  out  in  the  time  domain. 

The  numerical  Implementation  of  nonlinear  geometrical  acoustics 
is  divided  into  three  parts.  First  the  ray  paths  are  calculated  using  the  ray 
tracing  program  MEDUSA.  Next  the  program  CALC6  uses  the  environmental 
information  to  perform  a  numerical  intregration  along  the  ray  path, 
thereby  calculating  the  quantity  6.  The  ray  path  information,  the  value  6, 
and  the  environmental  Information  are  all  inputs  to  the  program  PLPROP. 
PLPROP  uses  this  information  to  calculate  the  distortion  range  variable,  Z. 
PLPROP  has  at  its  center  the  subroutine  WAVPROP,  a  finite  difference 
implementation  of  weak-shock  theory  for  plane  waves  (Pestorius  1973). 
Within  PLPROP  losses  are  accounted  for  in  the  frequency  domain  using  the 
empirical  relations  of  Frangois  and  Garrision  ( 1 982).  PLPROP  has  a  variety 
of  switches  that  allow  the  investigator  to  selectively  include  the  effects  of 
nonlinear  distortion,  ordinary  attenuation,  and  dispersion.  The  finite 
amplitude  effects  may  be  "turned  off"  after  a  specific  propagation  distance 
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has  been  reached.  The  output  of  the  program  is  in  the  form  of  plots  of  the 
time  waveforms  and  energy  spectra  for  preselected  propagation  distances. 
The  accuracy  of  the  program  was  verified  by  comparing  its  results  with 
those  of  known  analytic  solutions  for  waves  propagating  through  a  lossless 
homogeneous  media.  The  peak  pressure,  the  relative  shock  arrival  time,  and 
the  t/e  decay  time  of  a  weak  shock  with  an  exponentially  decaying  tail  were 
all  predicted  by  PLPROP  to  a  high  degree  of  accuracy 

In  the  numerical  study  the  effects  of  inhomogeneity,  ordinary 
attenuation  and  dispersion,  and  nonlinear  distortion  are  Investigated  by 
considering  the  propagation  of  explosion  waveforms  The  signals  used  in 
the  Investigation  are  a  weak  shock  with  an  exponentially  decaying  tail 
(exponential  pulse)  and  a  more  realistic  explosion  waveform  which  includes 
the  first  bubble  pulse.  The  exponential  pulse  is  propagated  58. 1  km  along  a 
ray  path  starting  at  a  depth  of  300  m.  The  more  realistic  waveform  is 
propagated  23  km  along  a  ray  path  starting  at  a  depth  of  4300  m. 

An  exponential  pulse  corresponding  to  a  0.818  kg  TNT  explosion  is 
used  to  Investigate  the  separate  and  combined  effects  of  inhomogeneity, 
ordinary  attenuation,  dispersion,  and  nonlinear  distortion.  The  effect  of 
Inhomogeneity  of  the  ocean  on  nonlinear  distortion  is  found  to  be  small. 
Dispersion  Is  found  to  play  an  Important  role  in  the  arrival  time  of  the  peak 
pressure.  The  propagation  of  a  finite  amplitude  signal  is  found  to  depend  on 
the  combination  of  nonlinear  distortion  and  ordinary  attenuation  and 
dispersion. 
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More  realistic  waveforms  corresponding  to  two  different 
explosions,  0.818  kg  and  22.7  kg  TNT,  are  used  to  quantify  the  answer  to  the 
Question,  "to  what  distance  are  finite  amplitude  effects  important?"  The 
spectra  of  the  signals  obtained  by  numerically  propagating  the  signals  to  a 
distance  of  23  km  accounting  for  ordinary  attenuation  and  dispersion  over 
the  entire  23  km  and  acounting  for  finite  amplitude  effects  as  indicated 
below: 

Case  A:  finite  amplitude  effects  neglected  entirely. 

Case  B:  finite  amplitude  effects  included  only  up  to  distance  150  m. 
Case  C;  finite  amplitude  effects  included  only  up  to  distance  1 1 00  m. 
Case  D;  finite  amplitude  effects  included  for  the  entire  23  km  path. 
The  Case  D  spectrum,  relative  to  that  of  Case  A,  is  lower  in  the  middle 
frequency  band  and  slightly  higher  in  the  low  and  high  frequency  bands.  A 
small  downward  shift  of  the  spectral  peaks  in  the  low  frequency  band  is 
also  noted.  For  the  one  ray  path  and  waveform  considered,  finite  amplitude 
effects  are  found,  in  the  case  of  a  0.818  kg  TNT  explosion,  to  be  of  small 
consequence  for  frequencies  below  6  kHz  and  distances  beyond  1 100  m.  For 
a  22.7  kg  TNT  explosion  the  corresponding  quantities  are  frequencies  below 
4  kHz  and  distances  beyond  1 1 00  m. 


APPENDIX  A 


RAY  COORDINATES  AND  THE  EXPRESSION  FOR  VH 


In  this  appendix  ray  coordinates  are  defined,  and  it  is  shown  that 
v2t=Ao'’  a(Ao/Co)/as,  where  Aq  is  the  area  of  the  ray  tube,  Cq  is  the 
small-signal  sound  speed,  and  s  is  the  distance  along  the  ray  path.  The 
expression  for  is  required  to  simplify  the  transport  equation, 

Eq.  (3-C.  1 9).  To  derive  the  expression,  we  use  tensor  analysis  and  a 
nonorthogonal  coordinate  system  called  ray  coordinates  (see,  for  example, 
Pitre  1984,  p.  54).  General  introductions  to  tensor  analysis  may  be  found  in 
standard  texts  such  as  those  by  Sokolnikoff  (1964)  and  Synge  and  Schild 
(1978). 

The  mathematical  development  proceeds  as  follows.  The  ray 
coordinate  system  and  the  Jacobian  of  the  transformation  from  Cartesian 
coordinates  to  ray  coordinates  are  found.  It  Is  then  shown  that  the 
determinant  of  the  Jacobian  is  related  to  the  ray  tube  area.  Next  the 
covariant  and  contravarlant  forms  of  the  metric  tensor  for  the  ray 
coordinate  system  are  found.  The  contravarlant  form  of  the  metric  tensor 
and  the  determinant  of  the  Jacobian  are  then  substituted  Into  the  definition 


of  the  Laplacian  for  a  general  coordinate  system.  The  desired  result  for 
follows  directly. 

We  define  some  notation:  Let  us  denote  the  determinant  of  a 
tensor  by  the  same  symbol,  but  without  indices;  for  example,  |A'i|  Is 
denoted  A.  Let  x‘,  where  I  =  l.  2, 3,  stand  for  the  rectangular  Cartesian 
coordinate  system,  that  Is,  x’  =  x,  x^  =  y,  x^-  z.  Let  stand  for  another  set 
of  coordinates,  ray  coordinates,  namely  x’  =  s,  =  <1>,  and  x^  =  where  s 
Is  the  distance  along  the  ray  path,  and  ^  and  are  the  ray  launch  angles 
with  respect  to  the  x  and  z  axes.  Both  the  ray  coordinate  system  and  the 
Cartesian  system  are  shown  In  Fig.  3.2.  Ray  coordinates  may  be  thought  of 
as  a  coordinate  system  based  on  the  path  of  a  particle  through  space 

A  physical  understanding  of  the  transformations  between  the  two 
coordinate  systems  Is  sought.  We  define  the  coordinate  transformations  as 
follows: 

x'  -  x'(x)  ,  (A.  1) 

and 

^  =  x''(x)  (A.2) 


The  Jacobian  JY  of  the  transformation  from  Cartesian  coordinates  to  ray 
coordinates  is 


i  _ 

'' "  ax'* 


dx/ds  dy/ds  dz/ds 
dx/d4>  dy/d(|)  dz/d4> 
dx/d<|>  dy/d<i>  dz/d>v. 


(A.3) 


Let  T  be  the  derivative  of  the  ray  path  position  vector  r  (see  Fig.  A.  1 ): 


FI6URE 


where  the  bracketed  terms  are  the  components  in  the  x,  y,  and  z  directions. 
The  derivatives  with  respect  to  the  launch  angles  <|>  and  denoted  and 
V^,  are  defined  as  follows; 


■  [d^’  d«l>’  d^l 


(A.5) 


\d«|»  dq>  d^f 


(A.6) 


Recall  from  vector  calculus  that  the  magnitude  of  the  vector  product, 

B  X  C,  equals  the  area  of  the  parallelogram  defined  by  B  and  C.  The  vector 
product  of  with  defines  a  cross-sectional  area  of  the  ray  tube  The 
scalar  triple  product,  T-(  x  V^),  defines  Aq,  the  cross-sectional  area  of 
the  ray  tube  normal  to  the  direction  of  propagation.  As  can  be  seen  from 
Eg.  (A.3),  the  scalar  triple  product  is  also  equal  to  the  determinant  of  the 
Jacobian; 


The  Laplactan  of  a  function  f  in  any  general  coordinate  system  is 
define J  as  follows  (Sokolnikoff  1964,  Eq.  92.1 1 ); 


where  f  is  a  function  in  the  original  coordinate  system.  To  use  Eq.  (A.8)  we 
need  the  contravariant  form  of  the  metric  tensor  in  the  ray  coordinate 
system.  To  find  we  first  find  the  covariant  form  Since  the 
contravariant  and  the  covariant  metric  tensors  are  inverses  of  one  another, 
the  contravariant  metric  tensor  can  he  found  from  the  covariant. 

In  a  Euclidian  space  of  three  dimensions,  the  covariant  forrn  of  the 
metric  tensor  may  be  defined  by  considering  an  elemental  distance  os  It  is 
known  (Sokolnikoff  1964,  Art.  29)  that  in  a  Cartesian  system,  the  square  of 
the  elemental  distance  ds  may  be  defined  as 

ds^  =  2  dx'  dx'  ,  (A.9) 

i 

or 

ds^  =  dx' dxj  ,  (A.  10) 

where  8jj  is  the  covariant  form  of  the  Kronecker  delta  8jj  for  the  Cartesian 
coordinate  system.  (In  the  Cartesian  coordinate  system,  the  contravariant 
and  covariant  forms  of  a  tensor  are  the  same,  therefore  8jj  =  8U  =  8jj ) 

We  seek  an  expression  for  dx'  to  use  in  Eq.  (A.9).  It  follows  from 


If  Eq.  (A.  1 1 )  Is  substituted  into  Eq.  (A.9),  it  can  be  seen  that 
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(A.  12) 


Equation  (A.  12)  may  be  written  as  follows  (Synge  and  Schild  1978,  Eq.  2. 102): 


where 


ds^  =  0^^  d/  dx*  . 


■?(#)(#) 


(A.  13) 


(A.  14) 


Equation  (A.  14)  Is  the  covariant  form  of  the  metric  tensor  for  the  ray 
coordinate  system. 

The  contravariant  form  of  the  metric  tensor  is  now  sought  The 
contravariant  and  covariant  forms  of  the  metric  tensor  are  inverses;  this 
may  be  expressed  mathematically  as  follows: 

9r,9"‘-S,‘  (A  15) 

Use  of  Eqs.  (A.  1 5)  and  (A.  1 4)  leads  to  the  following  expression  for  the 
contravariant  form  of  the  metric  tensor  (Guy  1985): 


•*I*S 

|SI . 


dx''\  /dx® 


We  need  to  know  the  determinant  of  the  metric  tensor.  From 
Eqs.  (A.3)  and  (A.  1 4),  it  can  be  shown  that 

-in 


(A.  1 7) 


The  general  form  of  the  Laplacian,  Eq.  (A.8),  can  now  be  calculated. 
Let  the  function  f  of  Eq.  (A.8)  be  the  eikonal  t.  Examination  of  the 
expression  (dt/dx^)  yields 


I.  .  / 

r"  [dx'jWl 


(A.  1 8] 


If  both  sides  are  multiplied  by  it  is  seen  that 


y  f  / dx*\ 

j  IdxV  IdxV  \dx7  WX 7 


(A  19) 


But  from  Eq.  (3-D.4)  it  is  known  that 


dxJ  Co\axVL 


It  therefore  follows  that 


U?  CoUxVUW 


=  Co-‘{1.0.0]  . 


(A21) 


If  Eqs.  (A.  16)  and  (A.2I )  are  substituted  into  Eq.  (A.8),  it  is  seen  that 


=  J-'  d(J/Co)/ds  . 


(A.22) 


Use  of  Eq.  (A.7)  yields 


VH  =  Aq-’  a(Ao/Co)/ds  , 


(A23) 


which  is  Eq.  (3-E.  1 ).  The  Laplacian  of  the  eikonal  is  proportional  to  the 
variation  of  the  ratio  of  ray  tube  area  to  the  small-signal  sound  speed. 
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APPENDIX  B 


THE  PARAMETER  OF  NONLINEARITY  FOR  SEAWATER 


In  this  appendix  the  analytic  formulation  of  the  parameter  of 
nonlinearity  B/A  (see,  for  example,  Beyer  1974,  p.  99)  Is  discussed  Some 
empirical  relations  developed  by  Morfey  (1984c)  for  use  In  the  numerical 
evaluation  of  the  parameter  for  seawater  are  then  presented.  The 
parameter  of  nonlinearity  accounts  for  the  curvature  of  the  pressure - 
density  relationship  of  a  fluid.  It  can  be  seen  from  Fig.  B.  1  that  by  moving  a 
finite  amount  from  the  quiescent  values  of  the  density  and  pressure,  and 
pQ,  the  slope  of  the  pressure-density  curve  changes.  By  definition,  the 
square  of  the  sound  speed  is  equal  to  the  slope  of  the  pressure-density 
curve, 


c 


2 


(B.l) 


An  Investigation  of  the  definition  of  the  sound  speed  yields  an  analytic 
expression  for  the  parameter  of  nonlinearity.  This  expression  may  be 
evaluated  numerically  by  use  of  empirical  formulas. 


In  order  to  develop  an  expression  for  c,  we  start  with  the 
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equation  of  state  for  seawater; 


P-P(p^X)  , 

where  %  is  the  salinity,  and  X  is  the  entropy.  Using  a  Taylor  series 
expansion  and  retaining  terms  up  to  second  order  yields 


(B.2) 


where 


x.t4)  -  pQ 


PoCo' 


B 


t^p  -^X.t4>*Po 


(B.3) 


(B.4) 


(B.5) 


n  =  <P  '  Po>/Po  ' 


(B.6) 


n  is  called  the  condensation,  in  arriving  at  Eg.  (B.3).  we  assumed  that 
variations  in  the  salinity  and  entropy  are  of  second  order  in  smallness. 
Equation  (B.  1 )  can  he  expanded  by  differentiating  Eq.  (B.3)  with  respect  to  p. 


r2  -  ?Dl 
■  l«n  »p]t.x 

»  (A  ♦  Bn)/pQ  ,  (B.7) 


»(A/Pq)(1  ♦qB/A)  .  (B.8) 


Using  a  binomial  expansion  and  retaining  terms  up  to  second  order  yields 
the  following  equation; 

c-(A/Po)''2(,  .  (B.9) 

The  next  step  is  to  determine  the  sound  speed  as  a  function  of  the 
particle  velocity  u.  in  the  case  of  progressive  plane  waves  traveling  in  the 
positive  X  direction  It  can  be  shown  that 

n*u/CQ  .  (B.10) 

The  use  of  Eg.  (B.  10)  In  the  case  of  outgoing  waves  In  an  inhomogeneous 
medium  deserves  some  attention.  In  the  derivation  of  nonlinear  acoustic 
theory  presented  In  Chapter  4,  we  used  ray  theory.  There  It  is  seen  that  an 
arbitrary  sound  wave  may  be  regarded  as  being  composed  of  many  locally 
plane  waves  (Landau  and  Lifshitz  1959,  p.  256).  Since  Eq.  (B.10)  is  valid  for 
plane  waves  and  the  wavefronts  in  the  inhomogeneous  medium  are  locally 
plane,  it  Is  assumed  that,  by  envoking  the  substitution  rule  (Chapter  2, 
Section  C),  Eq.  (B.10)  may  be  used  to  eliminate  n  trom  Eq.  (B.9).  Use  of 
Eq.  (B.10)  along  with  the  definition  of  A,  Eq.  (B.4),  yields  the  following 
equation; 

c  »  Cq  ♦  (B/2A)u  ,  (BID 

where  u  Is  the  particle  velocity  In  the  direction  of  the  ray  path.  In  Eq.  (B.  1 1 ) 
it  Is  seen  that  within  the  second-order  approximation,  the  sound  speed  of  a 
finite-amplitude  disturbance  Is  the  sum  of  the  small-signal  sound  speed  Cq 


and  a  perturbation  term  (B/2AXJ. 

We  now  recast  B/2A  in  a  way  that  lends  itself  to  numerical 
evaluation.  Use  of  the  definitions  of  A  and  B,  Eqs.  (B.4)  and  (B.5),  leads  to 


the  following: 


Pofrfr)] 


W4>“Po 


Algebraic  manipulation  of  terms  held  under  the  same  constraints  yields 


t.X4>*Po 


Following  the  method  of  Beyer  (1974,  p.  1 00),  we  use  the  chain  rule  to 
express  the  partial  derivative  of  Cq  in  terms  of  temperature  and  pressure. 


The  result  is  as  follows: 


^•P-P-P^ 


Where 


raT]  OtCT*  273.15) 
l^Pjx  P^P 


Qt  is  the  coefficient  of  thermal  expansion,  Cp  is  the  specific  heat  at 
constant  pressure,  and  T  is  the  temperature  in  ®C.  Morfey  (1984c)  has 


developed  empirical  relations  for  the  density  and  the  specific  heat  at 
constant  pressure.  These  relations  are  given  below.  Chen  and  Mi  Hero 
(1976)  have  developed  an  expression  for  Qj  which  is  valid  to  pressures  of 
100  jiPa.  Thus  B/2A  can  be  numerically  evaluated  over  a  wide  range  of 
salinities,  temperatures,  and  pressures. 

The  following  Is  Morfey’s  empirical  relation  for  the  density  of 
both  fresh  and  seawater,  p(P.T,U-  The  relative  error  is  typically  10'®'^  for 
the  temperature  range  0®  <  T  <  100®C  and  0.5  x  10"®^  for  the  range 
ICP  <  T  <  40®C.  The  Independent  variables,  their  units,  and  the  ranges  over 
which  they  may  vary  are  listed  below: 

P  »  pressure  (Pa) 

T  =  temperature  (®C) 

I  =  salinity  (»/^) 

The  empirical  equation  is  as  follows: 


p  -  1000/u  , 

(B.16) 

where 

u  =  Vq  ♦  ^/F  , 

(B.17) 

where 

(A,„.A„T)t  , 

(B.18) 

where 

F  =  5:,^,  3B,T'  ♦  B,o^  ♦P/(9.80665x  10^) 

+  6,(P/(9.80665x10^)]2  . 

(B.19) 

The  constants  are  defined  In  Table  B.  1. 


absolute  0<P<10^  , 

zero  salinity  0<T<100  , 

salt  water  0  <  T  <  40  , 

30  <  S  <  40  , 
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TABLE  B.  I 

CONSTANTS  REQUIRED  IN  THE  EMPIRICAL  EQUATION  FOR  DENSITY 


^  B,  Vi  6i 


0  0.25357  X 1 0*^  0.7 1 540  x  1 0*^  0.64575  X  1 0*^  -0.23430  x  1 0"^ 

1  0. 1 4894X  1 0*02  0.43 1 24  x  1 0*02 

2  -0.72023  X  1 0'O'  -0.36324  x  1 0*00 

3  -0.9 1030x1 0'O^  0.24620  x  I  O^o^ 

10  -0.35771x10*0'  0.62851  x10*0' 

11  -0.95 1 38  X  10^02 
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Morfey's  empirical  relation  for  the  specific  heat  at  constant  pressure 
Cp(ET.E)  is  given  below.  For  seawater  at  pressures  up  to  100  pPa,  we 
estimate  Cp  from 

Cp  *  C,  aX)  *  Cq  [P /(9.80665  x  10^)  +  HE,  T) 

-Cq(1+KE,T)  ,  (B.20) 

where  K  =  1 2  (bar/  is  an  empirical  constant  which  relates  the 
increase  of  Cp  with  pressure  for  seawater  to  that  for  freshwater.  The 
error  in  Cp  is  typically  lO"^"*.  The  function  C,  gives  the  specific  heat  of 
seawater  at  atmospheric  pressure  and  is  defined  as  follows; 


C,  =  Z ^  Z. B,.  T'’’  Ej’’ 

1  i=t.4  ij  ^ 


(B.21) 


The  error  in  this  equation  is  typically  lO"^  over  the  temperature  range 
0  -  40®  C  and  over  the  salinity  range  0-40  The  function  Cq  gives 
specific  heat  of  pure  water  and  is  defined  as  follows: 

'^0  “  ^1-1.4  ^j.1.3  l£^''/(9.80665x)0‘'))  (B.22) 

This  equation  is  valid  for  E  “  0  (zero  salinity)  over  the  temperature  and 
pressure  ranges  T  =  0  -  90®  C  and  P  -  0  -  100  pPa.  The  relative  error  is 
typically  10"®5.  The  constants  Ajj  and  Bjj  are  listed  in  Table  B.2. 
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TABLE  B.2 

CONSTANTS  REQUIRED  IN  THE  EMPIRICAL  EQUATION  FOR  SPECIFIC 
HEAT  AT  CONSTANT  PRESSURE 


i 

j 

Ai 

j 

BlJ 

1 

1 

0.42160 

XIO*®^ 

0.42179  XlO^O'* 

2 

1 

-0.24895 

xlO*®’ 

-0.34218 

3 

1 

0.46273 

xIO*^’ 

0.97816  xlO'^^’ 

4 

1 

-0.22446 

xIO'®^ 

-0.91609  x\0-^^ 

1 

2 

-0.46649 

xio+oo 

-0.72849  xlO*o> 

2 

2 

0.12204 

xIO'O’ 

0.19149  xlO*0® 

3 

2 

-0.19241 

xlO"®^ 

-0.65990  xl0-®2 

4 

2 

0.98686 

xIO"^ 

0.73197  xIO-^ 

1 

3 

0.12458 

xIO’^^^ 

0.22227  xlO-0' 

2 

3 

-0.25157 

X  10*^5 

-0.27630  X  10*02 

3 

3 

0.16518 

xIO'^^^ 

0.11588  X  10*05 

4 

3 

0.37555 

xlO’'° 

-0.13749  X  10*05 

APPENDIX  C 


ANALYTIC  SOLUTIONS  FOR  FINITE  AMPLITUDE  WAVES 
VIA  WEAK  SHOCK  THEORY 


In  this  appendix  the  analytic  solution  for  a  weak  shock  with  an 
exponentially  decaying  tail  propagating  through  a  homogeneous  medium  is 
developed.  The  solution  Is  used  to  verify  the  accuracy  of  the  propagation 
routine  discussed  in  Chapter  5;  the  results  of  the  numerical  routine  are 
compared  with  those  of  the  analytic  solution.  The  development  proceeds  as 
follows.  The  propagation  of  a  finite  amplitude  wave  through  a  homogeneous 
medium  is  considered.  A  general  solution  Is  obtained  for  the  case  In  which 
the  waveform  contains  weak  shocks.  This  solution  is  not  new;  see,  for 
example,  Blackstock  (1972).  Since  the  propagation  routine  operates  with 
nondlmenslonal  variables,  the  solution  is  nondlmensionalized.  The 
nondlmenslonal  solution  is  then  applied  to  the  particular  problem  of  the 


propagation  of  a  weak  shock  with  an  exponentially  decaying  tall.  The 
solution  for  this  waveform  Is  not  new  either  (Rogers  1977;  also  see 


Blackstock  1983). 


Before  discussing  the  propagation  of  weak  shocks,  we  examine  the 
propagation  of  a  continuous  (l.e.,  no  shocks  present)  finite  amplitude  wave 
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through  a  homogeneous  medium.  The  following  equation  is  the  exact  plane 
wave  equation  for  a  gas; 

u^  ♦  C(^Uj^  ♦  puu^  =  0  .  (C.  1 ) 

The  Earnshaw  solution  (see.  for  example,  Blackstock  1962,  Eq.  18)  satisfies 
Eq.  (C.  1 )  exactly.  However,  Eamshaw’s  solution  is  a  bit  cumbersome,  and 
particularly  so  for  problems  Involving  sources  for  which  the  boundary 
condition  is  given  in  the  form 

u  -  f(t)  at  X  *  0  .  (C.2) 

It  is  convenient  to  use  an  approximate  wave  equation.  Use  of  the 
first-order  plane  wave  equation  and  the  substitution  rule  of  Chapter  2 
converts  Eq.  (C.  I )  into 


U,  +  Cq\  -  pCp-^uu^  =  0 


(C.3) 


Use  of  the  definition  of  retarded  time  for  a  plane  wave,  Eq.  (3-C.l ), 
transforms  Eq.  (C.3)  into 


u^  -  PCq-^uu^.  -  0 


(C.4) 


An  equally  valid  expression  In  terms  of  the  presssure  P'  is 

PV(pW)P’PV  -0  • 


(C.5) 


Equation  (C.5)  Is  obtained  from  Eq.  (C.4)  by  using  the  expression  (see,  for 
example,  Blackstock  1962,  Eq.  95) 


A 


U  =  P'/PflCo  -  pPV2P(,V 


(C.6) 


The  "Earnshaw  solution"  of  Eq.  (C.6)  for  the  boundary  condition 

P' =  g(t)  =  g(t')  at  x  =  0 
is 

P'  =  gi^) 

where 

♦  =  t' ♦  pg(4*)x/poCo3 


(C.7) 

(C.8) 

(C.Q) 


By  using  the  transforms  of  the  independent  and  dependent  variables  given  in 
Table  C.  I,  the  solution  presented  in  Eqs.  (C.8)  and  (C.9)  can  be  used  for 
nonplanar  geometries  and  inhomogeneous  fluids  (Blackstock  1966,  Carlton 
and  Blackstock  1974). 

TABLE  C.l 

DISTORTION  RANGE  AND  TRANSFORM  PRESSURE  VARIABLES 


Plane 

Cylindrical 

Spherical 

Inhomogeneous  Fluid 


Z 

X 

Cq  An  (r/r^) 
Sq  in  (Gs/Sq) 


W 

P' 

P’  (r/r^)’'^ 

P'  (r/r^j) 

P'  <ApoCo/Aopc)’'2 


Here  r  is  the  range  and  ro  is  the  reference  range.  If  the  boundary  condition 


W  -  g(t)  at  Z  -  0 


(CIO) 


the  general  solution  to  a  second  approximation  is  as  follows; 

W-g(^)  .  (C.11) 

where 

♦  -  t‘ ♦  g(^)pZ/poCo3  .  (Cl 2) 

This  solution  is  referred  to  as  the  generalized  approximate  Eamshaw 
solution. 

We  now  discuss  weak  shocks  within  the  context  of  the  Eamshaw 
solution.  The  approximate  Eamshaw  solution  is  valid  on  either  side  of  a 
weak  shock;  knowledge  of  the  shock  strength  allows  one  to  tie  the  two 
Eamshaw  solutions  together.  The  relative  arrival  time  of  the  shock  is, 
however,  unknown;  we  therefore  seek  an  equation  for  the  relative  shock 
arrival  time,  it  turns  out,  however,  that  we  must  first  find  an  expression 
for  the  shock  velocity. 

Application  of  the  equations  of  hydrodynamics  to  a  propagating 
weak  shock  leads  to  an  equation  for  the  shock  velocity.  Consider  a  plane 
weak  shock  propagating  through  a  homogeneous  fluid.  The  fluid  ahead  of  the 
shock  is  assumed  to  be  undisturbed;  its  properties  are  denoted  by  the 
subscript  a.  The  properties  of  the  fluid  behind  the  shock  are  denoted  by  the 
subscript  b.  If  the  continuity,  momentum,  and  energy  equations  are  applied 
to  the  control  volume  moving  with  the  shock  velocity  V^,  several  relations, 
called  the  Rankine-Hugoniot  shock  relations,  can  be  derived: 


PbfV*  -  =  P.v*  . 

V  -  PS^sh  -  V  “  0 


Following  Blackstock  ( 1 972),  we  solve  Eq.  (C.  1 5)  for  an  approximate 
expression  for  the  shock  velocity 

Vsh  =  Co*PV2 


(C  13) 


(C.14) 


(C.15) 


(C.16) 


If  the  particle  velocity  in  front  of  the  wave  had  been  the  nonzero  value  u,, 
the  following  equation  would  have  obtained; 


Vsh  *  S  * 


(C.17) 


Since  the  particle  velocity  both  ahead  of  and  behind  the  shock  are  known 
from  the  two  Earnshaw  solutions,  the  shock  velocity  can  be  calculated. 

Use  of  the  expression  for  the  shock  velocity,  Eq.  (C.  17),  leads  to  an 
expression  for  the  relative  shock  arrival  time.  Still  following  Blackstock, 
we  assume  that  the  shock  first  occurred  at  time  1  and  position  n.  The 
current  position  and  time  of  the  shock  are  denoted  by  Xjj,  and  tg^, 
respectively,  where 


(C  18) 


If  a  binomial  expansion  is  used  and  terms  up  to  second  order  are  retained, 
Eq.  (C.  1 7)  may  be  expressed  as 


Substitution  of  Eq.  (C.  19)  Into  Eq.  (CIS)  and  simplification  of  the  result 
using  the  definition  of  the  retarded  time  variable,  Eq.  (3-C.  1 ),  as  well  as  the 
progressive  wave  Impedance  relation,  yields 

f*  - 1'  -  (P/2q^Cq%  (p;  ♦  PV  dx  .  (C.20) 

The  above  equation  can  be  expressed  in  the  terms  of  the  transformed 
independent  and  dependent  variables,  W  and  Z,  since  the  Ranklne-Hugonlot 
relations  are  invariant  under  the  tranformatlons  (Blackstock  1966, 

Appendix  A).  Restating  Eq.  (C.20)  In  terms  of  W  and  Z  gives 

(p/2PoCo’)Ij(W,«W,)dZ  .  (C.2I) 

The  differential  form  of  Ea  (C.21)  is 

dt'3^/dZ  •  -p(W^  ♦  W,)/2poCo3  .  (C.22) 

Equations  (C.21)  and  (C.22)  give  the  relative  shock  arrival  time.  Other 
useful  relations  regarding  the  relative  shock  arrival  time  are  obtained 
from  the  solution  for  the  propagation  of  continuous  waves,  Eqs.  (C.  1 1 )  and 
(C.12).  By  noting  that  the  value  of  W  is  different  before  and  after  the  shock, 
we  obtain  the  following ; 


l'*  ■  ♦b  -  (PZ/Po^o’*  V  >  23) 


t'rt  ■  ♦.  -  (pZ/Po^o’)  9(*.) 


(C.24) 


As  was  mentioned  earlier,  we  are  developing  this  solution  to 
compare  with  the  results  of  the  propagation  routine.  Since  the  computer 
algorithm  operates  on  nondimensional  variables,  we  nondimensionalize  the 
solution.  This  is  done  using  characteristic  dimensions  of  the  problem.  The 
nondimensional  variables  are  defined  below: 


V  =  W/P‘o  . 

(C.25) 

t  a  t/tc  , 

(C.26) 

r'atVt,  . 

(C.27) 

0  S  pWoZ/PoCo^t^  , 

(C.28) 

§  =  «/tc  , 

(C.29) 

eaP'/PoCo^  . 

(C.30) 

G($) «  g(^)/P’o  , 

(C.31) 

where  t,.  is  a  characteristic  time,  and  P'q  is  the  reference  pressure 
amplitude.  In  terms  of  the  nondimensional  variables,  the  nonlinear  wave 
equation  and  its  solution  are  expressed  as  follows. 


Nonlinear  wave  equation 

< 

o 

1 

< 

< 

II 

O 

(C.32) 

Boundary  condition 

V  =  G(t)  at  0  *  0 

(C.33) 

Continuous  wave  solution 

V  =  G(§)  . 

(C.34) 

where 

$  ■  t'  ♦  0G(5)  , 

(C.35) 

Relative  shock  arrival  time 

(C.36) 

We  now  use  the  nondimenslonal  solution  to  solve  for  the  relative 


shock  arrival  time,  the  shock  amplitude,  and  the  I/e  decay  time  of  a  weak 
shock  with  an  exponentially  decaying  tail.  The  boundary  condition  for  the 


waveform  is  as  follows: 


which  converts  to 


P'  =  A  exp  (-t/lQ)  H(t)  at  r  =  r^  , 


W  «  A  exp  (-t/T.)  H(t)  at  Z  -  0  , 


(C.30) 


(C.39) 


where  H(t)  is  the  step  function.  Expressing  the  boundary  condition  in  terms 
of  the  nondimenslonal  variables  yields 

V  *  exp  (-t)  H(t)  at  0  *  0  ,  (C.40) 

where  the  characteristic  values  t^  and  Wq  are  defined  to  be  and  A, 
respectively.  Use  of  Eqs.  (C.34)  and  (C.35)  leads  to  the  solution  to  the 
continuous  portion  of  wave: 


where 


V  =  exp  (-$)  , 

$  *  r’  ♦  0  exp  (-$) 


(C.41) 


(C.42) 


I 

4 


We  now  solve  for  the  shock  amplitude.  Considering  only  the  shock, 
we  may  write  Eqs  (C.41 )  and  (C.42)  as  follows: 


Taking  the  natural  logarithm  of  Eq.  (C.43),  and  then  substituting  the  result 
into  Eq.  (C.44)  and  rearranging,  we  obtain  the  following: 

.  (C.45) 

By  taking  the  derivative  of  Eq.  (C.45)  with  respect  to  o  and  using  Eq.  (C.37) 
noting  that  V3  is  equal  to  zero,  we  arrive  at  the  following  expression; 

-1/2  *  (1  ♦  0V^)  dV^/do  “  0  .  (C.46) 

If  Eq.  (C.46)  is  integrated  with  respect  to  0,  noting  that  the  integration 
constant  must  be  chosen  to  fit  the  boundary  conditions,  a  quadratic  equation 
with  the  following  roots  is  obtained: 

V^  =  [-l  i  7(1  *20)1/0  .  (C.47) 

The  positive  root  in  Eq.  (C.47)  allows  the  boundary  condition,  V,,  =  I  at  a  =0, 

to  be  recovered. 

We  now  solve  for  the  relative  shock  arrival  time  We  start  by 
substituting  the  positive  root  of  Eq.  (C.47)  into  Eq.  (C  36),  again  noting  that 
Vg  is  zero.  If  the  starting  time  of  the  shock  is  assumed  to  be  0,  the  result 
is  as  follows: 

t’s^  = -UnV^  *  oVbl  .  (C48) 

The  expression  for  the  1/e  decay  time  is  obtained  using  Eqs  (r  47) 


and  (C  48).  The  1/e  decay  time  is  defined  as  the  time  required  for  the  pulse 
to  decay  to  1/e  of  its  current  value,  V|,/e.  Substitution  of  V^/e  for  Vj,  in 
Eg.  (C.48)  yields  the  relative  time  of  the  1/e  point  of  the  waveform.  The 
following  expression  for  the  1/e  decay  time  is  obtained  by  subtracting  the 
relative  time  of  1/e  point  from  the  relative  shock  arrival  time, 

i:',^,  =  -[l +0  Vj,(l  -  1/e)]  .  (C.49) 

where  r'  is  the  1  /e  decay  time. 

In  summary  we  have  found  the  solution  for  the  shock  amplitude, 

the  relative  shock  arrival  time,  and  the  1/e  decay  time  of  a  weak  shock 
with  an  exponentially  decaying  tail,  Eqs.  (C.47),  (C.48),  and  (C.49), 
respectively.  This  information,  coupled  with  the  definitions  of  W  and  Z 
described  in  Table  C.l,  permits  us  to  calculate  the  shock  amplitude,  relative 
shock  arrival  time,  and  the  1/e  decay  time  independent  of  whether  the  wave 
is  plane,  cylindrical,  spherical,  or  confined  to  a  ray  tube  in  an 
inhomogeneous  medium. 


APPENDIX  D 


COMPUTER  PROORAM 

Computer  program  PLPROP  is  a  FORTRAN  4  program  which 
implements  nonlinear  geometrical  acoustics.  The  program  is  described  at 
length  in  Chapter  5.  The  heart  of  the  program  is  the  pair  of  subroutines 
WAVPROP  and  RESAMP  which  were  written  by  Pestorius  ( 1 973).  The  fast 
Fourier  Transform  routine  and  plotting  routines,  which  were  written  at 
Applied  Research  Laboratories,  The  University  of  Texas  at  Austin,  Austin, 
Texas,  are  used  extensively  in  the  program. 

The  input  files  and  other  parameters  required  of  the  user  are 
described  in  the  program.  The  computer  code  has  a  running  commentary  on 
its  function  which  should  help  the  user  understand  operation  of  the 
program. 
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